library(R.matlab)
library(rlang)
library(languageR)
library(lme4)
library(MASS)
library(ggplot2) #cookbook for R/ Graphs
library(hexbin)
library(memisc)
library(reshape)
library(reshape2) #melt and cast -> restructure and aggregate data
library(data.table)
library(coin) #for permutation tests
library(psych)
library(doBy)
library(heplots)
library(plyr) #necessary for ddply
library(matrixStats) 
library(foreign) 
library(Hmisc)
#library(lmerTest) # fucks with tab_model in sjPlot, so I don't use that anymore.
library (stringr)
library(gridExtra)
library(grid)
library(gdata)
library(effects)
library(ggExtra)
library("piecewiseSEM")
library(ordinal)
library(simr)
library(sjPlot)
library(htmlTable)
library(lmerTest)
detach("package:lmerTest", unload = TRUE)
setwd("/Users/Romy/Dropbox (Brown)/ShenhavLab/Experiments/bas/Analysis/HDDM/LMMGLMR")
source("/Users/Romy/Beruflich/R-funcs/summarySEwithinO.R")
source("/Users/Romy/Beruflich/R-funcs/summarySE.R")
source("/Users/Romy/Beruflich/R-funcs/normDataWithin.R")
source("/Users/Romy/Beruflich/R-funcs/multiplot.R")

Loading all data

### Study 1: choose best/ choose worst behavioral
input_file = '/Users/Romy/Dropbox (Brown)/ShenhavLab/Experiments/bas/Analysis/HDDM/BASB_108Data_041118.xls'
a1 = read.xls(input_file)


### Study 2: fMRI
input_file = '~/Dropbox (Brown)/ShenhavLab/Experiments/bas/Analysis/BASB_508Data_040918.xls' # 508 data

a1fMRI = read.xls(input_file)

## BARTRA ROI
input_file = '~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/BASBWfMRI/BASB_508BARTRAROIData_053118.xls'

BAROI = read.xls(input_file)

## vStr and vmPFC from BARTRA
input_file = '~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/BASBWfMRI/BASB_508BARTRAROIsubRegData_112118.xls'

BAROIsub = read.xls(input_file)

## SEED ROIs within vStr
input_file = '~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/BASBWfMRI/BASB_508_left_righ_Str_ROIS.xls'

SEEDROIS = read.xls(input_file)

## ROI of OVr activation
input_file = '~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/BASBWfMRI/BASB_508_OVr_ROI.xls'

OVrROI = read.xls(input_file)

# Shenhav & Buckner (2014) ROIs for PCC 
input_file = '~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/BASBWfMRI/BASB_508ROIData_053118.xls' 

ROI = read.xls(input_file)

# vmPFC subregions Shenhav & Karmarkar (2019)
input_file = '~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/BASBWfMRI/VMPFC_OFC_ROI.xls' 

VMPFCROI = read.xls(input_file)

### Study 0: choose best only
input_file = '/Users/Romy/Dropbox (Brown)/ShenhavLab/Experiments/bas/Analysis/HDDM/BASB_105Data_150318.xls' 
a0 = read.xls(input_file)

### Study 3: choose worst only
input_file = '/Users/Romy/Dropbox (Brown)/ShenhavLab/Experiments/bas/Analysis/HDDM/BASB_106Data_03222018.xls' 
a2 = read.xls(input_file)

### Study 4: choose best/ choose worst - incentivized
input_file = '/Users/Romy/Dropbox (Brown)/ShenhavLab/Experiments/bas/Analysis/BASB_118Data_072619.xls'

a118 = read.xls(input_file)

Study 1: best vs worst within subject

data prep

# 1: center covariates
a1$cgoalvAvRemBid <- scale(a1$goalvAvRemBid/10, center=TRUE, scale = FALSE)                  ## unsigned value difference (max/min vs average remaining)
a1$cgoalvNext <- scale(a1$goalvNext, center=TRUE, scale = FALSE)
a1$cAvBid <- scale(a1$AvBid/10, center=TRUE, scale = FALSE)
a1$cRT <- scale(a1$logRT, center=TRUE, scale=FALSE)
#a1$cChosenvAvRemV <- scale(a1$ChosenvAvRemV, center= TRUE, scale=FALSE)
a1$Liking <- scale(a1$tmpAllLrateAllTrials, scale=FALSE, center=TRUE)
a1$avGoalValue <- a1$AvBid                                                                ## the value of the set relative to the goal
a1$avGoalValue[a1$isChooseBestTrial==0] <- (a1$AvBid [a1$isChooseBestTrial==0]*-1) +10
a1$cavGoalValue <- scale(a1$avGoalValue/10, scale=FALSE, center= TRUE)
a1$goalChosenV <- a1$ChosenV                                                              ## the value of the chosen Item relativ to goal
a1$goalChosenV[a1$isChooseBestTrial==0] <- (a1$ChosenV[a1$isChooseBestTrial==0]*-1) +10
a1$cgoalChosenV <- scale(a1$goalChosenV, scale=FALSE, center= TRUE)
a1$cChosenV <- scale(a1$ChosenV, scale=FALSE, center=TRUE)
a1$avUnchV <- a1$SumUnchV/3                                                               ## average of the unchosen values
a1$cavUnchV <- scale(a1$avUnchV, scale=FALSE, center=TRUE)
a1$avgoalUnchV <- a1$avUnchV                                                              ## average of unchosen values relative to goal
a1$avgoalUnchV[a1$isChooseBestTrial==0] <- (a1$avUnchV[a1$isChooseBestTrial==0]*-1) +10
a1$cavgoalUnchV <- scale(a1$avgoalUnchV, scale=FALSE, center=TRUE)
a1$relGoalvAvRem <- a1$goalChosenV - a1$avgoalUnchV                                       ## goal related chosen vs unchosen (negative Values are suboptimal choices) 
a1$crelGoalvAvRem <- scale(a1$relGoalvAvRem/10, scale=FALSE, center=TRUE)
a1$goalVal <- a1$MaxBid
a1$goalVal[a1$isChooseBestTrial==0] <- a1$MinBid[a1$isChooseBestTrial==0]
a1$cgoalVal <- scale(a1$goalVal, scale=FALSE, center=TRUE)                               ## goal Value (best matching value in set relative to goal)
a1$rewSpgoalvAvRem <- a1$goalVal - a1$avUnchV                                            ## goal v av rem in bid space
a1$crewSpgoalvAvRem <- scale(a1$rewSpgoalvAvRem/10, scale=FALSE, center=TRUE) 
a1$Prgoalcons <- a1$goalVal
a1$Prgoalcons[a1$isChooseBestTrial==0] <- 10 - a1$MinBid[a1$isChooseBestTrial==0]        ## this is a terrible variable name, but that's the revalued goal value
a1$cChosenvAvRemV <- scale(a1$ChosenvAvRemV/10, scale=FALSE, center=TRUE)   
a1$chosenvsnext <- a1$ChosenV - a1$allTrProdBid2

# 2: make categorical predictors factors and set contrasts
a1$subj_idx <- as.factor(a1$subj_idx)
a1$isChooseBestTrial <- as.factor(a1$isChooseBestTrial)
contrasts(a1$isChooseBestTrial) <- contr.sdif(2)
a1$block <- rep(1, length(a1$subj_idx))
a1$block[a1$Trial > 60] <- 2
a1$block <- as.factor(a1$block)
contrasts(a1$block) <- contr.sdif(2)
a1$CondType <- as.factor(a1$CondType)
a1$CondType <- revalue(a1$CondType, c("1"="mixed", "2"="low", "3"="medium", "4"="high"))


## exclude trials with RT longer than 15s
as <- a1
a1 <- a1[a1$rt<15,]

accuracy

We analyzed choice consistency as a function of Overall Value, Value Difference, Task goal, and also tested for an interaction of Value Difference with Task goal, to see whether participants were more sensitive to value information in either condition.

BWaccmod0t1 <- glmer(response~ cgoalvAvRemBid*isChooseBestTrial+cAvBid+(isChooseBestTrial|subj_idx), data = a1[!a1$goalvNext==0,], #
                                  family = binomial) 

sv1_max <- svd(getME(BWaccmod0t1, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(BWaccmod0t1, transform= NULL ,show.stat = TRUE, string.stat="z",   dv.labels = c("P(Goal Chosen)"), pred.labels= c("(Intercept)", "Value Difference","Best - Worst Conditon", "Overall Value", "Value Difference:Best - Worst Conditon"), col.order = c("est", "se","ci",  "stat","df", "p"))

	P(Goal Chosen)
Predictors	Log-Odds	CI	z	p
(Intercept)	-0.22	-0.37 – -0.07	-2.91	0.004
Value Difference	3.76	3.20 – 4.32	13.12	<0.001
Best - Worst Conditon	0.11	-0.19 – 0.40	0.69	0.488
Overall Value	0.33	-0.27 – 0.92	1.08	0.280
Value Difference:Best - Worst Conditon	0.58	-0.54 – 1.70	1.01	0.311
Random Effects
σ²	3.29
τ₀₀ _{subj_idx}	0.03
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.13
ρ₀₁ _{subj_idx}	-0.07
ICC	0.02
N _{subj_idx}	30
Observations	1703
Marginal R² / Conditional R²	0.169 / 0.183

There was no significant interaction, so we excluded the interaction term.

BWaccmod0 <- glmer(response~ cgoalvAvRemBid+isChooseBestTrial+cAvBid+(isChooseBestTrial|subj_idx), data = a1[!a1$goalvNext==0,], #
                                  family = binomial)

sv1_max <- svd(getME(BWaccmod0, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(BWaccmod0, transform= NULL ,show.stat = TRUE, string.stat="z",   dv.labels = c("P(Goal Chosen)"), pred.labels= c("(Intercept)", "Value Difference","Best - Worst Conditon", "Overall Value"), col.order = c("est", "se","ci",  "stat","df", "p"))

	P(Goal Chosen)
Predictors	Log-Odds	CI	z	p
(Intercept)	-0.22	-0.37 – -0.07	-2.90	0.004
Value Difference	3.75	3.19 – 4.30	13.15	<0.001
Best - Worst Conditon	0.19	-0.06 – 0.44	1.52	0.129
Overall Value	0.29	-0.29 – 0.88	0.98	0.327
Random Effects
σ²	3.29
τ₀₀ _{subj_idx}	0.03
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.13
ρ₀₁ _{subj_idx}	-0.24
ICC	0.02
N _{subj_idx}	30
Observations	1703
Marginal R² / Conditional R²	0.167 / 0.181

eff_df <- Effect(c("cgoalvAvRemBid", "isChooseBestTrial"),BWaccmod0, xlevels=list(cgoalvAvRemBid=seq(min(a1$cgoalvAvRemBid), max(a1$cgoalvAvRemBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))
IA$cgoalvAvRemBid <- (IA$cgoalvAvRemBid - min(IA$cgoalvAvRemBid))*10

## with CI instead of se
pBWacc <- ggplot(data=IA, aes(x=cgoalvAvRemBid, y=fit , color=isChooseBestTrial )) + geom_line()+scale_colour_manual(name="Goal",values=c("#CC0000","#000000")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cgoalvAvRemBid, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous( expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("value difference") + ylab("Pr(goal chosen)") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.justification=c(1,0), legend.position=c(0.95,0.05))
pBWacc

pdf("/Users/Romy/Dropbox (Brown)/ShenhavLab/Experiments/bas/Analysis/HDDM/best_worst_within_main_accuracy_w_choice_condition_as_color.pdf", width = 4, height = 4)#, units = 'cm', res = 200, compression = 'lzw'
pBWacc 
dev.off()

pdf 
  2

S1BWaccmod0G <- glmer(response~ cgoalvAvRemBid*isChooseBestTrial+cavGoalValue+(isChooseBestTrial|subj_idx), data = a1[!a1$goalvNext==0,], 
                                  family = binomial) 

sv1_max <- svd(getME(S1BWaccmod0G, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(S1BWaccmod0G, transform= NULL ,show.stat = TRUE, string.stat="z",   dv.labels = c("P(Goal Chosen)"), pred.labels= c("(Intercept)", "Value Difference","Best - Worst Conditon", "Overall Goal Value", "Value Difference:Best - Worst Conditon"), col.order = c("est", "se","ci",  "stat","df", "p"))

	P(Goal Chosen)
Predictors	Log-Odds	CI	z	p
(Intercept)	-0.21	-0.36 – -0.06	-2.77	0.006
Value Difference	3.72	3.16 – 4.29	12.92	<0.001
Best - Worst Conditon	0.11	-0.19 – 0.41	0.74	0.459
Overall Goal Value	-0.27	-0.87 – 0.32	-0.90	0.369
Value Difference:Best - Worst Conditon	0.50	-0.62 – 1.62	0.87	0.386
Random Effects
σ²	3.29
τ₀₀ _{subj_idx}	0.03
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.13
ρ₀₁ _{subj_idx}	0.00
ICC	0.02
N _{subj_idx}	30
Observations	1703
Marginal R² / Conditional R²	0.170 / 0.186

In summary, choice consistency was only sensitive to value difference.
***

rt

Basic checks

# make table of RT by task

meanRTbySubS1 <- ddply(a1, .(subj_idx, isChooseBestTrial), summarise,
                  mrt    = median(rt, na.rm=TRUE))

meanRTDiff <-  meanRTbySubS1$mrt[meanRTbySubS1$isChooseBestTrial==0] -  meanRTbySubS1$mrt[meanRTbySubS1$isChooseBestTrial==1]

ggplot(meanRTbySubS1, aes(isChooseBestTrial,mrt,  color= isChooseBestTrial)) + geom_violin()+geom_boxplot() + theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))

ggplot(a1, aes(rt, group =isChooseBestTrial, color=isChooseBestTrial, fill=isChooseBestTrial)) + geom_density(alpha=0.5) + theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))

run analyses separately for best and worst before putting all in one model

BWRTmod1bBEST <- lmer(logRT~ cgoalvAvRemBid+cAvBid+(cgoalvAvRemBid|subj_idx), data = a1[a1$isChooseBestTrial==1,],
                                     REML=FALSE)

sv1_max <- svd(getME(BWRTmod1bBEST, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  


BWRTmod1bWORST <- lmer(logRT~ cgoalvAvRemBid+cAvBid+(cgoalvAvRemBid|subj_idx), data = a1[a1$isChooseBestTrial==0,],
                                      REML=FALSE)

sv1_max <- svd(getME(BWRTmod1bWORST, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(BWRTmod1bBEST,BWRTmod1bWORST, p.val= "kr", show.df = T, transform= NULL ,show.stat = TRUE, string.stat="t",   dv.labels = c("log RT Best", "log RT worst"),pred.labels= c("(Intercept)", "Value Difference", "Overall Value"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT Best					log RT worst
Predictors	Estimates	CI	t	df	p	Estimates	CI	t	df	p
(Intercept)	1.39	1.31 – 1.46	36.34	31.00	<0.001	1.46	1.38 – 1.53	38.01	31.00	<0.001
Value Difference	-0.50	-0.60 – -0.39	-9.37	29.00	<0.001	-0.49	-0.59 – -0.40	-10.02	30.00	<0.001
Overall Value	-0.33	-0.42 – -0.25	-7.72	1719.00	<0.001	0.42	0.34 – 0.50	10.39	1714.00	<0.001
Random Effects
σ²	0.18					0.16
τ₀₀	0.04 _{subj_idx}					0.04 _{subj_idx}
τ₁₁	0.02 _{subj_idx.cgoalvAvRemBid}					0.02 _{subj_idx.cgoalvAvRemBid}
ρ₀₁	-0.57 _{subj_idx}					-0.03 _{subj_idx}
ICC	0.18					0.20
N	30 _{subj_idx}					30 _{subj_idx}
Observations	1737					1733
Marginal R² / Conditional R²	0.084 / 0.249					0.102 / 0.280

All data: standard analysis with Overall and Relative Value, plus main effect, but no interaction with choice goal

# all data no interaction
BWRTmod0 <- lmer(logRT~ cgoalvAvRemBid+cAvBid+isChooseBestTrial+(isChooseBestTrial|subj_idx), data = a1,
                                 REML=FALSE)

sv1_max <- svd(getME(BWRTmod0, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(BWRTmod0,p.val= "kr", show.df = T, transform= NULL ,show.stat = TRUE, string.stat="t",   dv.labels = c("log RT"),pred.labels= c("(Intercept)", "Value Difference", "Overall Value", "Best - Worst Condition"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.42	1.36 – 1.49	42.11	31.00	<0.001
Value Difference	-0.46	-0.52 – -0.40	-15.19	3432.00	<0.001
Overall Value	0.05	-0.01 – 0.11	1.65	3447.00	0.099
Best - Worst Condition	-0.06	-0.13 – 0.01	-1.81	31.00	0.081
Random Effects
σ²	0.18
τ₀₀ _{subj_idx}	0.03
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.03
ρ₀₁ _{subj_idx}	0.01
ICC	0.18
N _{subj_idx}	30
Observations	3470
Marginal R² / Conditional R²	0.058 / 0.224

All data: adding interaction of Overall Value with Choice goal

# all data with interaction
BWRTmod1b <- lmer(logRT~ cgoalvAvRemBid+cAvBid*isChooseBestTrial+(isChooseBestTrial|subj_idx), data = a1,
                                 REML=FALSE)

sv1_max <- svd(getME(BWRTmod1b, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(BWRTmod1b,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("log RT"), pred.labels= c("(Intercept)", "Value Difference", "Overall Value", "Best - Worst Condition", "Overall Value: Best - Worst"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.42	1.36 – 1.49	41.93	31.00	<0.001
Value Difference	-0.49	-0.55 – -0.43	-16.48	3430.00	<0.001
Overall Value	0.05	-0.01 – 0.11	1.61	3448.00	0.108
Best - Worst Condition	-0.07	-0.14 – 0.00	-1.90	31.00	0.066
Overall Value: Best - Worst	-0.75	-0.86 – -0.63	-12.66	3296.00	<0.001
Random Effects
σ²	0.17
τ₀₀ _{subj_idx}	0.03
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.03
ρ₀₁ _{subj_idx}	-0.00
ICC	0.18
N _{subj_idx}	30
Observations	3470
Marginal R² / Conditional R²	0.096 / 0.262

Nested model to get Overall Value effects for each choice goal

## nested model
BWRTmod1bn <- lmer(logRT~ isChooseBestTrial/(cAvBid)+cgoalvAvRemBid+(isChooseBestTrial|subj_idx), data = a1,
                                  REML=FALSE)

sv1_max <- svd(getME(BWRTmod1bn, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  


tab_model(BWRTmod1bn,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("log RT"), pred.labels= c("(Intercept)", "Best - Worst Condition", "Value Difference", "Worst: Overall Value", "Best: Overall Value"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.42	1.36 – 1.49	41.93	31.00	<0.001
Best - Worst Condition	-0.07	-0.14 – 0.00	-1.90	31.00	0.066
Value Difference	-0.49	-0.55 – -0.43	-16.48	3430.00	<0.001
Worst: Overall Value	0.42	0.34 – 0.50	10.09	3418.00	<0.001
Best: Overall Value	-0.33	-0.41 – -0.24	-7.76	3391.00	<0.001
Random Effects
σ²	0.17
τ₀₀ _{subj_idx}	0.03
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.03
ρ₀₁ _{subj_idx}	-0.00
ICC	0.18
N _{subj_idx}	30
Observations	3470
Marginal R² / Conditional R²	0.096 / 0.262

Figure 1 C overlay part 1

eff_df <- Effect(c("cAvBid", "isChooseBestTrial"), BWRTmod1b)
IA <- as.data.frame(eff_df)
IA$cAvBid <- (IA$cAvBid - min(IA$cAvBid))*10
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
pASVBW <- ggplot(data=IA, aes(x=cAvBid, y=fit , color=isChooseBestTrial )) + geom_line()+scale_colour_manual(name="Goal",values=c("#CC0000","#000000")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cAvBid, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous(limits=c(0.8, 2.1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Overall value") + ylab(" ") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+ theme(legend.justification=c(1,1), legend.position=c(0.95,0.25) )


eff_df <- Effect(c("cgoalvAvRemBid", "isChooseBestTrial"), BWRTmod1b, xlevels=list(cgoalvAvRemBid=seq(min(a1$cgoalvAvRemBid), max(a1$cgoalvAvRemBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))
IA$cgoalvAvRemBid <- (IA$cgoalvAvRemBid - min(IA$cgoalvAvRemBid))*10
## with CI instead of se
pVDBW <- ggplot(data=IA, aes(x=cgoalvAvRemBid, y=fit , color=isChooseBestTrial )) + geom_line()+scale_colour_manual(name="Goal",values=c("#CC0000","#000000")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cgoalvAvRemBid, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous(limits=c(0.8, 2.1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Relative Value") + ylab("log RT") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+ theme(legend.position="none" )
multiplot(pVDBW, pASVBW,  cols=2)

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/HDDM/best_worst_within_allRS.pdf", width = 8, height = 4)#, units = 'cm', res = 200, compression = 'lzw'
multiplot(pVDBW, pASVBW,  cols=2)
dev.off()

pdf 
  2

Modeling goal congruency instead of reward value

## goal value instead of AvBid
BWRTmod1g <- lmer(logRT~ cgoalvAvRemBid+cavGoalValue+isChooseBestTrial+(isChooseBestTrial|subj_idx), data = a1,
                                 REML=FALSE)

sv1_max <- svd(getME(BWRTmod1g, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  


tab_model(BWRTmod1g,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("log RT"), pred.labels= c("(Intercept)", "Value Difference", "Overall Goal Value", "Best - Worst Condition"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.42	1.36 – 1.49	42.01	31.00	<0.001
Value Difference	-0.49	-0.55 – -0.43	-16.42	3429.00	<0.001
Overall Goal Value	-0.37	-0.43 – -0.32	-12.66	3298.00	<0.001
Best - Worst Condition	-0.08	-0.15 – -0.01	-2.21	31.00	0.035
Random Effects
σ²	0.17
τ₀₀ _{subj_idx}	0.03
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.03
ρ₀₁ _{subj_idx}	0.00
ICC	0.18
N _{subj_idx}	30
Observations	3470
Marginal R² / Conditional R²	0.095 / 0.261

eff_df <- Effect(c("cavGoalValue", "isChooseBestTrial"), BWRTmod1g)
IA <- as.data.frame(eff_df)
IA$cavGoalValue <- (IA$cavGoalValue - min(IA$cavGoalValue))*10
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
pASVBW <- ggplot(data=IA, aes(x=cavGoalValue, y=fit , color=isChooseBestTrial )) + geom_line()+scale_colour_manual(name="Goal",values=c("#CC0000","#000000")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cavGoalValue, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous(limits=c(0.8, 2.1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Overall value") + ylab(" ") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.justification=c(1,1), legend.position=c(0.95,0.25))

eff_df <- Effect(c("cgoalvAvRemBid", "isChooseBestTrial"), BWRTmod1g, xlevels=list(cgoalvAvRemBid=seq(min(a1$cgoalvAvRemBid), max(a1$cgoalvAvRemBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$cgoalvAvRemBid <- (IA$cgoalvAvRemBid - min(IA$cgoalvAvRemBid))*10
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
pVDBW <- ggplot(data=IA, aes(x=cgoalvAvRemBid, y=fit , color=isChooseBestTrial )) + geom_line()+scale_colour_manual(name="Goal",values=c("#CC0000","#000000")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cgoalvAvRemBid, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous(limits=c(0.8, 2.1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Relative Value") + ylab("log RT") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="none" )
multiplot( pVDBW,pASVBW, cols=2)

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/HDDM/best_worst_within_allGS.pdf", width = 8, height = 4)#, units = 'cm', res = 200, compression = 'lzw'
multiplot( pVDBW,pASVBW, cols=2)
dev.off()

pdf 
  2

test for residual OV effects

BWRTmod1allin <- lmer(logRT~ cgoalvAvRemBid+cavGoalValue+ cAvBid+isChooseBestTrial+(isChooseBestTrial|subj_idx), data = a1,
                                 REML=FALSE)

sv1_max <- svd(getME(BWRTmod1allin, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  


tab_model(BWRTmod1allin,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("log RT"), pred.labels= c("(Intercept)", "Value Difference", "Overall Goal Value", "Overall Reward Value", "Best - Worst Condition"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.42	1.36 – 1.49	41.97	31.00	<0.001
Value Difference	-0.49	-0.55 – -0.43	-16.48	3430.00	<0.001
Overall Goal Value	-0.37	-0.43 – -0.32	-12.66	3296.00	<0.001
Overall Reward Value	0.05	-0.01 – 0.11	1.61	3448.00	0.108
Best - Worst Condition	-0.08	-0.15 – -0.01	-2.21	31.00	0.035
Random Effects
σ²	0.17
τ₀₀ _{subj_idx}	0.03
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.03
ρ₀₁ _{subj_idx}	-0.00
ICC	0.18
N _{subj_idx}	30
Observations	3470
Marginal R² / Conditional R²	0.096 / 0.262

Model comparison

Table 2 part 1

# simplest model with only VD
BWRTmodbase <- lmer(logRT~ cgoalvAvRemBid+isChooseBestTrial+(isChooseBestTrial|subj_idx), data = a1,
                                     REML=FALSE)

sv1_max <- svd(getME(BWRTmodbase, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  


tab_model(BWRTmodbase,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("log RT"), pred.labels= c("(Intercept)", "Value Difference", "Best - Worst Condition"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.42	1.36 – 1.49	42.16	31.00	<0.001
Value Difference	-0.46	-0.52 – -0.40	-15.13	3431.00	<0.001
Best - Worst Condition	-0.06	-0.13 – 0.01	-1.81	31.00	0.080
Random Effects
σ²	0.18
τ₀₀ _{subj_idx}	0.03
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.03
ρ₀₁ _{subj_idx}	0.01
ICC	0.18
N _{subj_idx}	30
Observations	3470
Marginal R² / Conditional R²	0.058 / 0.223

xxx <- anova(BWRTmod0, BWRTmod1b, BWRTmod1g, BWRTmodbase)
row.names(xxx) <- c("VD + C<sub>b-w</sub> (baseline)", "baseline + OV<sub>reward</sub>", "baseline + OV<sub>goal</sub>", "baseline + C<sub>b-w</sub>: OV<sub>reward</sub>")
attributes(xxx)$heading <- NULL
xxx2 <- xxx#[, c(1,2, 6, 7, 8)]
xxx2$AIC <- round(xxx2$AIC)
xxx2$BIC <- round(xxx2$BIC)
xxx2$Chisq <- round(xxx2$Chisq,2)
xxx2[,8] <- round(xxx2[,8],3)
xxx2$logLik <- xxx2$AIC - c(NaN, xxx2$AIC[1:length(xxx2$AIC)-1]) 
allR2s <- rsquared(list(BWRTmodbase,BWRTmod0, BWRTmod1b, BWRTmod1g))
xxx2[,1] <- round(allR2s$Conditional,2)
xxx2 <- xxx2[, c(1:4, 6, 8)]
colnames(xxx2) <- c("R<sup>2</sub>", "AIC", "BIC", "dAIC", "Chi<sup>2</sub>", "p")
htmlTable(xxx2)

	R²	AIC	BIC	dAIC	Chi²	p
VD + C_b-w (baseline)	0.22	4112	4156
baseline + OV_reward	0.23	4112	4161	0	2.72	0.099
baseline + OV_goal	0.26	3957	4006	-155	155.09	0
baseline + C_b-w: OV_reward	0.26	3956	4011	-1	2.59	0.108

Liking

Liking is related to reward values, not goal values

a1$LikingOrdinal <- as.factor(a1$tmpAllLrateAllTrials)
#ordered(a1fMRI$LikingOrdinal, levels=c("1", "2", "3", "4", "5"))
a1$cMaxV <- scale(a1$MaxBid/10, scale=FALSE, center=TRUE)
fmm1Lik <- clmm(LikingOrdinal ~ isChooseBestTrial+crelGoalvAvRem +cAvBid +     cavGoalValue +cRT+(isChooseBestTrial+cRT|subj_idx), data= a1[!a1$Choice1==-1,])

tab_model(fmm1Lik, transform= NULL,show.stat = TRUE, string.stat = "z", pred.labels= c("(Intercept: 1|2)","(Intercept: 2|3)", "(Intercept: 3|4)", "(Intercept: 4|5)","Best - Worst Condition", "Value Difference","Overall Reward Value", "Overall Goal Value",  "RT"))

	LikingOrdinal
Predictors	Log-Odds	CI	z	p
(Intercept: 1\|2)	-2.67	-2.97 – -2.37	-17.26	<0.001
(Intercept: 2\|3)	-0.52	-0.81 – -0.23	-3.54	<0.001
(Intercept: 3\|4)	1.30	1.01 – 1.59	8.78	<0.001
(Intercept: 4\|5)	3.16	2.85 – 3.47	20.05	<0.001
Best - Worst Condition	0.07	-0.11 – 0.25	0.73	0.465
Value Difference	-0.11	-0.36 – 0.13	-0.92	0.357
Overall Reward Value	6.92	6.57 – 7.26	39.09	<0.001
Overall Goal Value	-0.10	-0.37 – 0.17	-0.73	0.467
RT	-0.18	-0.36 – -0.00	-1.98	0.048
N _{subj_idx}	30
Observations	3469

fmm1Likb <- clmm(LikingOrdinal ~ isChooseBestTrial+cMaxV +cAvBid +     cavGoalValue +cRT+(isChooseBestTrial+cRT|subj_idx), data= a1[!a1$Choice1==-1,])

tab_model(fmm1Likb, transform= NULL,show.stat = TRUE, string.stat = "z", pred.labels= c("(Intercept: 1|2)","(Intercept: 2|3)", "(Intercept: 3|4)", "(Intercept: 4|5)","Best - Worst Condition", "Max Value","Overall Reward Value", "Overall Goal Value",  "RT"))

	LikingOrdinal
Predictors	Log-Odds	CI	z	p
(Intercept: 1\|2)	-2.69	-2.99 – -2.38	-17.21	<0.001
(Intercept: 2\|3)	-0.51	-0.80 – -0.22	-3.44	0.001
(Intercept: 3\|4)	1.31	1.02 – 1.61	8.78	<0.001
(Intercept: 4\|5)	3.15	2.84 – 3.46	19.89	<0.001
Best - Worst Condition	0.07	-0.11 – 0.25	0.73	0.463
Max Value	0.61	0.24 – 0.99	3.19	0.001
Overall Reward Value	6.35	5.86 – 6.84	25.44	<0.001
Overall Goal Value	-0.07	-0.35 – 0.20	-0.53	0.599
RT	-0.11	-0.29 – 0.07	-1.19	0.234
N _{subj_idx}	30
Observations	3469

## interlude: power analysis for Study 4 (incentivized version) We need 3-ish participants for 80% power. We chose to get 14 to be safe.
# power <- powerSim(BWRTmod1b,test = fixed("cAvBid:isChooseBestTrial"), nsim = 200)
# pc1 <- powerCurve(BWRTmod1b,test = fixed("cAvBid:isChooseBestTrial"), nsim = 200, along= "subj_idx")
# plot(pc1)

LCA simulations:

We simulated data from two LCA models using the behavioral data from Study 1. One took the original ratings as inputs (reward values) and the other one goal congruency coded values (goal values). Only the goal value LCA could produce the observed reversing overall value effect on RT.

Goal value LCA

input_file = '/Users/Romy/Dropbox (Brown)/ShenhavLab/Experiments/bas/Analysis/HDDM/allSubDataTableSIMBASBWflip_01152019.xls' # LCA simulated data

a4= read.xls(input_file)

a4$subj_idx <- as.factor(a4$subj_idx)
a4$cgoalvAvRemBid <- scale(a4$goalvAvRemBid, center=TRUE, scale = FALSE)
a4$cAvBid <- scale(a4$AvBid, center=TRUE, scale = FALSE)
a4$avGoalValue <- a4$AvBid  ## the value of the set relative to the goal
a4$avGoalValue[a4$isChooseBestTrial==0] <- (a4$AvBid[a4$isChooseBestTrial==0]*-1) +10
a4$cavGoalValue <- scale(a4$avGoalValue, scale=FALSE, center= TRUE)
a4$isChooseBestTrial <- as.factor(a4$isChooseBestTrial)
contrasts(a4$isChooseBestTrial) <- contr.sdif(2)

lSIMRT <- log(a4$SIMRT)
BWKRTmod0worstSIM <- lm(lSIMRT~ isChooseBestTrial*(cgoalvAvRemBid+cAvBid), data = a4)

Figure 1D

eff_df <- Effect(c("cAvBid", "isChooseBestTrial"), BWKRTmod0worstSIM,xlevels=list(cAvBid=seq(min(a4$cAvBid), max(a4$cAvBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$cAvBid <- IA$cAvBid +5
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
pgOV <- ggplot(data=IA, aes(x=cAvBid, y=fit , color=isChooseBestTrial )) + geom_line(size=2)+scale_colour_manual(name="Goal",values=c("#CC0000","#999999")) +theme_bw(12)+ 
  scale_y_continuous(limits=c(0.8, 1.65), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Overall Value") + ylab("log RT") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")

eff_df <- Effect(c("cgoalvAvRemBid", "isChooseBestTrial"), BWKRTmod0worstSIM,xlevels=list(cgoalvAvRemBid=seq(min(a4$cgoalvAvRemBid), max(a4$cgoalvAvRemBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$cgoalvAvRemBid <- IA$cgoalvAvRemBid - min(a4$cgoalvAvRemBid)
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
pgRV <- ggplot(data=IA, aes(x=cgoalvAvRemBid, y=fit , color=isChooseBestTrial )) + geom_line(size=2)+scale_colour_manual(name="Goal",values=c("#CC0000","#999999")) +theme_bw(12)+ 
  scale_y_continuous(limits=c(0.8, 1.65), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Value Difference") + ylab("log RT") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")

multiplot( pgOV,pgRV, cols=2)

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/HDDM/SIMbw_flip.pdf", width = 8, height = 4)
multiplot( pgOV,pgRV, cols=2)
dev.off()

pdf 
  2

SI Figures 2-3

## goal value model
BWKRTmod0worstSIM <- lm(lSIMRT~ isChooseBestTrial+(cgoalvAvRemBid+cavGoalValue), data = a4)


eff_df <- Effect(c("cavGoalValue", "isChooseBestTrial"), BWKRTmod0worstSIM,xlevels=list(cAvBid=seq(min(a4$cAvBid), max(a4$cAvBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$cavGoalValue <- IA$cavGoalValue +5
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
pgOV <- ggplot(data=IA, aes(x=cavGoalValue, y=fit , color=isChooseBestTrial )) + geom_line(size=2)+scale_colour_manual(name="Goal",values=c("#CC0000","#999999")) +theme_bw(12)+ 
  scale_y_continuous(limits=c(0.8, 1.65), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Overall Goal Value") + ylab("log RT") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")

pgOV

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/HDDM/SIMbw_GVRT.pdf", width = 4, height = 4)
pgOV
dev.off()

pdf 
  2

eff_df <- Effect(c("cgoalvAvRemBid", "isChooseBestTrial"), BWKRTmod0worstSIM,xlevels=list(cgoalvAvRemBid=seq(min(a4$cgoalvAvRemBid), max(a4$cgoalvAvRemBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$cgoalvAvRemBid <- IA$cgoalvAvRemBid - min(a4$cgoalvAvRemBid)
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
pgRV <- ggplot(data=IA, aes(x=cgoalvAvRemBid, y=fit , color=isChooseBestTrial )) + geom_line(size=2)+scale_colour_manual(name="Goal",values=c("#CC0000","#999999")) +theme_bw(12)+ 
  scale_y_continuous(limits=c(0.8, 1.65), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Value Difference") + ylab("log RT") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")

eff_df <- Effect(c("isChooseBestTrial"), BWKRTmod0worstSIM)
Catmain <- as.data.frame(eff_df)
Catmain$isChooseBestTrial<- revalue(Catmain$isChooseBestTrial, c("0"="worst", "1"="best"))
pGVmain <- ggplot(data=Catmain, aes(x=isChooseBestTrial, y=fit, color = isChooseBestTrial))+ geom_point(size=8, shape = 18)+scale_y_continuous(limits=c(0.8, 1.65), expand=c(0, 0))+
  geom_errorbar(aes(max = upper, min = lower), width=0,alpha=0.3) + geom_point()+ theme_bw()+theme(axis.text = element_text(size = 12), strip.text.x= element_text(size=12))+xlab("Choice Condition") +scale_colour_manual(name="Goal",values=c("#CC0000","#999999")) + ylab("log RT")+
  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")

multiplot( pgRV,pGVmain, cols=2)

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/HDDM/SIMbw_GV.pdf", width = 8, height = 4)
multiplot( pgRV,pGVmain, cols=2)
dev.off()

pdf 
  2

## plot for GV and ACC
BWKRTmod0worstSIM <- lm(SIMCHOICE~ isChooseBestTrial*(cgoalvAvRemBid+cavGoalValue), data = a4[!a4$goalvNext==0,])

eff_df <- Effect(c("cavGoalValue", "isChooseBestTrial"), BWKRTmod0worstSIM,xlevels=list(cAvBid=seq(min(a4$cAvBid), max(a4$cAvBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$cavGoalValue <- IA$cavGoalValue +5
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

pgCorrgv <- ggplot(data=IA, aes(x=cavGoalValue, y=fit , color=isChooseBestTrial )) + geom_line(size=2)+scale_colour_manual(name="Goal",values=c("#CC0000","#999999")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cavGoalValue, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous(limits=c(0, 1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Overall Value") + ylab("P(Goal Chosen)") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom") + ggtitle("Goal Value LCA")

eff_df <- Effect(c("cgoalvAvRemBid", "isChooseBestTrial"), BWKRTmod0worstSIM,xlevels=list(cgoalvAvRemBid=seq(min(a4$cgoalvAvRemBid), max(a4$cgoalvAvRemBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$cgoalvAvRemBid <- IA$cgoalvAvRemBid - min(a4$cgoalvAvRemBid)
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
pgRVACCgv <- ggplot(data=IA, aes(x=cgoalvAvRemBid, y=fit , color=isChooseBestTrial )) + geom_line(size=2)+scale_colour_manual(name="Goal",values=c("#CC0000","#999999")) +theme_bw(12)+ 
  scale_y_continuous(limits=c(0, 1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Value Difference") + ylab("P(Goal Chosen)") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")

Reward value LCA

SI Figures 1

input_file = '/Users/Romy/Dropbox (Brown)/ShenhavLab/Experiments/bas/Analysis/HDDM/allSubDataTableSIMBASBWreward_01152019.xls' # LCA simulated data

a4= read.xls(input_file)

a4$subj_idx <- as.factor(a4$subj_idx)
a4$cgoalvAvRemBid <- scale(a4$goalvAvRemBid, center=TRUE, scale = FALSE)
a4$cAvBid <- scale(a4$AvBid, center=TRUE, scale = FALSE)
a4$avGoalValue <- a4$AvBid  ## the value of the set relative to the goal
a4$avGoalValue[a4$isChooseBestTrial==0] <- (a4$AvBid[a4$isChooseBestTrial==0]*-1) +10
a4$cavGoalValue <- scale(a4$avGoalValue, scale=FALSE, center= TRUE)
a4$isChooseBestTrial <- as.factor(a4$isChooseBestTrial)
contrasts(a4$isChooseBestTrial) <- contr.sdif(2)

lSIMRT <- log(a4$SIMRT)
BWKRTmod0worstSIM <- lm(lSIMRT~ isChooseBestTrial*(cgoalvAvRemBid+cAvBid), data = a4)


eff_df <- Effect(c("cAvBid", "isChooseBestTrial"), BWKRTmod0worstSIM,xlevels=list(cAvBid=seq(min(a4$cAvBid), max(a4$cAvBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$cAvBid <- IA$cAvBid +5
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
prewOV <- ggplot(data=IA, aes(x=cAvBid, y=fit , color=isChooseBestTrial )) + geom_line(size=2)+scale_colour_manual(name="Goal",values=c("#CC0000","#999999")) +theme_bw(12)+ 
  scale_y_continuous(limits=c(0.8, 1.65), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Overall Value") + ylab("log RT") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")

eff_df <- Effect(c("cgoalvAvRemBid", "isChooseBestTrial"), BWKRTmod0worstSIM,xlevels=list(cgoalvAvRemBid=seq(min(a4$cgoalvAvRemBid), max(a4$cgoalvAvRemBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$cgoalvAvRemBid <- IA$cgoalvAvRemBid - min(a4$cgoalvAvRemBid)
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
prewRV <- ggplot(data=IA, aes(x=cgoalvAvRemBid, y=fit , color=isChooseBestTrial )) + geom_line(size=2)+scale_colour_manual(name="Goal",values=c("#CC0000","#999999")) +theme_bw(12)+ 
  scale_y_continuous(limits=c(0.8, 1.65), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Value Difference") + ylab("log RT") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")

multiplot(prewOV, prewRV, cols=2)

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/HDDM/SIMbw_rew.pdf", width = 8, height = 4)
multiplot(prewOV, prewRV, cols=2)
dev.off()

pdf 
  2

BWKRTmod0worstSIM <- lm(SIMCHOICE~ isChooseBestTrial*(cgoalvAvRemBid+cavGoalValue), data = a4[!a4$goalvNext==0,])

eff_df <- Effect(c("cavGoalValue", "isChooseBestTrial"), BWKRTmod0worstSIM,xlevels=list(cAvBid=seq(min(a4$cAvBid), max(a4$cAvBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$cavGoalValue <- IA$cavGoalValue +5
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

pgCorrrv <- ggplot(data=IA, aes(x=cavGoalValue, y=fit , color=isChooseBestTrial )) + geom_line(size=2)+scale_colour_manual(name="Goal",values=c("#CC0000","#999999")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cavGoalValue, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous(limits=c(0, 1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Overall Value") + ylab("P(Goal Chosen)") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")+ ggtitle("Reward Value LCA")

eff_df <- Effect(c("cgoalvAvRemBid", "isChooseBestTrial"), BWKRTmod0worstSIM,xlevels=list(cgoalvAvRemBid=seq(min(a4$cgoalvAvRemBid), max(a4$cgoalvAvRemBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$cgoalvAvRemBid <- IA$cgoalvAvRemBid - min(a4$cgoalvAvRemBid)
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
pgRVACCrv <- ggplot(data=IA, aes(x=cgoalvAvRemBid, y=fit , color=isChooseBestTrial )) + geom_line(size=2)+scale_colour_manual(name="Goal",values=c("#CC0000","#999999")) +theme_bw(12)+ 
  scale_y_continuous(limits=c(0, 1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Value Difference") + ylab("P(Goal Chosen)") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")

Accuracy of choices as a function of relative and overall goal value for the goal vs reward value LCA

multiplot(pgCorrgv, pgRVACCgv, pgCorrrv, pgRVACCrv, cols=2)

Study 2: best vs worst within fMRI

data prep

# 1: center covariates
a1fMRI$cgoalvAvRemBid <- scale(a1fMRI$goalvAvRemBid/10, center=TRUE, scale = FALSE)                  ## unsigned value difference (max/min vs average remaining)
a1fMRI$cgoalvNext <- scale(a1fMRI$goalvNext/10, center=TRUE, scale = FALSE)
a1fMRI$cAvBid <- scale(a1fMRI$AvBid/10, center=TRUE, scale = FALSE)
a1fMRI$cRT <- scale(a1fMRI$logRT, center=TRUE, scale=FALSE)
#a1fMRI$cChosenvAvRemV <- scale(a1fMRI$ChosenvAvRemV, center= TRUE, scale=FALSE)
a1fMRI$Liking <- a1fMRI$tmpAllLrateAllTrials
a1fMRI$avGoalValue <- a1fMRI$AvBid                                                               ## the value of the set relative to the goal
a1fMRI$avGoalValue[a1fMRI$isChooseBestTrial==0] <- (a1fMRI$AvBid [a1fMRI$isChooseBestTrial==0]*-1) +10
a1fMRI$cavGoalValue <- scale(a1fMRI$avGoalValue/10, scale=FALSE, center= TRUE)
a1fMRI$goalChosenV <- a1fMRI$ChosenV                                                              ## the value of the chosen Item relativ to goal
a1fMRI$goalChosenV[a1fMRI$isChooseBestTrial==0] <- (a1fMRI$ChosenV[a1fMRI$isChooseBestTrial==0]*-1) +10
a1fMRI$cgoalChosenV <- scale(a1fMRI$goalChosenV/10, scale=FALSE, center= TRUE)
a1fMRI$cChosenV <- scale(a1fMRI$ChosenV/10, scale=FALSE, center=TRUE)
a1fMRI$avUnchV <- a1fMRI$SumUnchV/3                                                               ## average of the unchosen values
a1fMRI$cavUnchV <- scale(a1fMRI$avUnchV/10, scale=FALSE, center=TRUE)
a1fMRI$avgoalUnchV <- a1fMRI$avUnchV                                                              ## average of unchosen values relative to goal
a1fMRI$avgoalUnchV[a1fMRI$isChooseBestTrial==0] <- (a1fMRI$avUnchV[a1fMRI$isChooseBestTrial==0]*-1) +10
a1fMRI$cavgoalUnchV <- scale(a1fMRI$avgoalUnchV/10, scale=FALSE, center=TRUE)
a1fMRI$relGoalvAvRem <- a1fMRI$goalChosenV - a1fMRI$avgoalUnchV                                       ## goal related chosen vs unchosen (negative Values are suboptimal choices) 
a1fMRI$crelGoalvAvRem <- scale(a1fMRI$relGoalvAvRem/10, scale=FALSE, center=TRUE)
a1fMRI$goalVal <- a1fMRI$MaxBid
a1fMRI$goalVal[a1fMRI$isChooseBestTrial==0] <- a1fMRI$MinBid[a1fMRI$isChooseBestTrial==0]
a1fMRI$cgoalVal <- scale(a1fMRI$goalVal, scale=FALSE, center=TRUE)                               ## goal Value (best matching value in set relative to goal)
a1fMRI$rewSpgoalvAvRem <- a1fMRI$goalVal - a1fMRI$avUnchV                                            ## goal v av rem in bid space
a1fMRI$crewSpgoalvAvRem <- scale(a1fMRI$rewSpgoalvAvRem/10, scale=FALSE, center=TRUE) 
a1fMRI$Prgoalcons <- a1fMRI$goalVal
a1fMRI$Prgoalcons[a1fMRI$isChooseBestTrial==0] <- 10 - a1fMRI$MinBid[a1fMRI$isChooseBestTrial==0]        ## this is a terrible variable name, but that's the revalued goal value
a1fMRI$ChosenvAvRemV <- scale(a1fMRI$ChosenvAvRemV/10, scale=FALSE, center=TRUE)    
a1fMRI$cLiking <- scale(a1fMRI$Liking/5, scale=FALSE, center=TRUE)

## max and min V in reward space
a1fMRI$cMaxVR <- scale(a1fMRI$MaxBid/10, scale=FALSE, center=TRUE)
a1fMRI$cMinVR <- scale(a1fMRI$MinBid/10, scale=FALSE, center=TRUE)

a1fMRI$MaxVG <- a1fMRI$MaxBid
a1fMRI$MaxVG [a1fMRI$isChooseBestTrial==0] <- 10 - a1fMRI$MaxBid[a1fMRI$isChooseBestTrial==0]
a1fMRI$MinVG <- a1fMRI$MinBid 
a1fMRI$MinVG [a1fMRI$isChooseBestTrial==0] <- 10 - a1fMRI$MinBid[a1fMRI$isChooseBestTrial==0]

a1fMRI$cMaxVG <- scale(a1fMRI$MaxVG/10, scale=FALSE, center=TRUE)
a1fMRI$cMinVG <- scale(a1fMRI$MinVG/10, scale=FALSE, center=TRUE)

# 2: make categorical predictors factors and set contrasts
a1fMRI$subj_idx <- as.factor(a1fMRI$subj_idx)
a1fMRI$isChooseBestTrial <- as.factor(a1fMRI$isChooseBestTrial)
contrasts(a1fMRI$isChooseBestTrial) <- contr.sdif(2)

a1fMRI$block <- rep(1, length(a1fMRI$subj_idx))
a1fMRI$block[a1fMRI$Trial > 60] <- 2
a1fMRI$block <- as.factor(a1fMRI$block)
contrasts(a1fMRI$block) <- contr.sdif(2)
a1fMRI$CondType <- as.factor(a1fMRI$CondType)
a1fMRI$CondType <- revalue(a1fMRI$CondType, c("1"="mixed", "2"="low", "3"="medium", "4"="high"))

## add and transform neural ROI data

#colnames(BAROI)
a1fMRI$binConjunc_PvNxDECxRECxMONxPRI <- scale(BAROI$binConjunc_PvNxDECxRECxMONxPRI,center = TRUE, scale= TRUE)

## with scaling
a1fMRI$sBARTRA_NAcc <- scale(BAROIsub$BARTRA_NAcc,center = TRUE, scale= TRUE)
a1fMRI$sBARTRA_vmPFC <- scale(BAROIsub$BARTRA_vmPF,center = TRUE, scale= TRUE)

a1fMRI$sDMN_Str <- scale(SEEDROIS$DMN_Str_4mm,center = TRUE, scale= TRUE)
a1fMRI$sFCN_Str <- scale(SEEDROIS$FCN_Str_4mm,center = TRUE, scale= TRUE)
a1fMRI$sLN_Str <- scale(SEEDROIS$LN_Str_4mm,center = TRUE, scale= TRUE)

a1fMRI$sPOS_PCC <- scale(ROI$POS_PCC_p001_k50_isolated__0__31_41,center = TRUE, scale= TRUE)

# BA14 [mOFC] vs. BA24 [pgACC])
a1fMRI$sHUvC_mOFC <- scale(VMPFCROI$Anat_ProbAtlas_BA14_p50bin_1,center = TRUE, scale= TRUE)
a1fMRI$sPOS_rACC <- scale(VMPFCROI$Anat_ProbAtlas_BA24_p50bin_1,center = TRUE, scale= TRUE)
## without scaling

a1fMRI$BARTRA_NAcc <- BAROIsub$BARTRA_NAcc
a1fMRI$BARTRA_vmPFC <- BAROIsub$BARTRA_vmPF

a1fMRI$POS_PCC <- ROI$POS_PCC_p001_k50_isolated__0__31_41
# BA14 [mOFC] vs. BA24 [pgACC])
a1fMRI$HUvC_mOFC <-VMPFCROI$Anat_ProbAtlas_BA14_p50bin_1
a1fMRI$POS_rACC <- VMPFCROI$Anat_ProbAtlas_BA24_p50bin_1

a1fMRI$DMN_Str <- SEEDROIS$DMN_Str_4mm
a1fMRI$FCN_Str <- SEEDROIS$FCN_Str_4mm
a1fMRI$LN_Str <- SEEDROIS$LN_Str_4mm

a1fMRI$OVrROI <- OVrROI$OVr_L_R
## exclude trials with RT longer than 15s
as <- a1fMRI
a1fMRI <- a1fMRI[a1fMRI$rt<15,]

accuracy

S2BWaccmod0 <- glmer(response~ cgoalvAvRemBid*isChooseBestTrial+cAvBid+(isChooseBestTrial|subj_idx), data = a1fMRI[!a1fMRI$goalvNext==0,], 
                                  family = binomial) 

sv1_max <- svd(getME(S2BWaccmod0, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(S2BWaccmod0, transform= NULL ,show.stat = TRUE, string.stat="z",   dv.labels = c("P(Goal Chosen)"), pred.labels= c("(Intercept)", "Value Difference","Best - Worst Conditon", "Overall Value", "Value Difference:Best - Worst Conditon"), col.order = c("est", "se","ci",  "stat","df", "p"))

	P(Goal Chosen)
Predictors	Log-Odds	CI	z	p
(Intercept)	-0.28	-0.49 – -0.07	-2.59	0.010
Value Difference	3.65	3.12 – 4.18	13.40	<0.001
Best - Worst Conditon	-0.01	-0.27 – 0.24	-0.10	0.918
Overall Value	0.49	-0.08 – 1.06	1.70	0.090
Value Difference:Best - Worst Conditon	0.53	-0.53 – 1.58	0.98	0.326
Random Effects
σ²	3.29
τ₀₀ _{subj_idx}	0.23
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.05
ρ₀₁ _{subj_idx}	-0.76
ICC	0.07
N _{subj_idx}	30
Observations	2136
Marginal R² / Conditional R²	0.144 / 0.202

eff_df <- Effect(c("cgoalvAvRemBid", "isChooseBestTrial"),S2BWaccmod0, xlevels=list(cgoalvAvRemBid=seq(min(a1fMRI$cgoalvAvRemBid), max(a1fMRI$cgoalvAvRemBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))
IA$cgoalvAvRemBid <- IA$cgoalvAvRemBid - min(IA$cgoalvAvRemBid)

## with CI instead of se
pBWaccBWcol <- ggplot(data=IA, aes(x=cgoalvAvRemBid, y=fit , color=isChooseBestTrial )) + geom_line()+scale_colour_manual(name="Goal",values=c("#CC0000","#000000")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cgoalvAvRemBid, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous( limits=c(0, 1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 1), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("value difference") + ylab("Pr(goal chosen)") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.justification=c(1,0), legend.position=c(0.95,0.05))

pBWaccBWcol

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/BASBWfMRI/best_worst_within_main_accuracyBWcolor.pdf", width = 4, height = 4)#, units = 'cm', res = 200, compression = 'lzw'
pBWaccBWcol 
dev.off()

pdf 
  2

S2BWaccmod0G <- glmer(response~ cgoalvAvRemBid*isChooseBestTrial+cavGoalValue+(isChooseBestTrial|subj_idx), data = a1fMRI[!a1fMRI$goalvNext==0,], 
                                  family = binomial) 

sv1_max <- svd(getME(S2BWaccmod0G, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(S2BWaccmod0G, transform= NULL ,show.stat = TRUE, string.stat="z",   dv.labels = c("P(Goal Chosen)"), pred.labels= c("(Intercept)", "Value Difference","Best - Worst Conditon", "Overall Goal Value", "Value Difference:Best - Worst Conditon"), col.order = c("est", "se","ci",  "stat","df", "p"))

	P(Goal Chosen)
Predictors	Log-Odds	CI	z	p
(Intercept)	-0.27	-0.48 – -0.07	-2.58	0.010
Value Difference	3.66	3.12 – 4.19	13.38	<0.001
Best - Worst Conditon	-0.00	-0.26 – 0.25	-0.02	0.982
Overall Goal Value	0.13	-0.39 – 0.64	0.48	0.632
Value Difference:Best - Worst Conditon	0.46	-0.59 – 1.50	0.85	0.395
Random Effects
σ²	3.29
τ₀₀ _{subj_idx}	0.22
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.04
ρ₀₁ _{subj_idx}	-0.71
ICC	0.06
N _{subj_idx}	30
Observations	2136
Marginal R² / Conditional R²	0.140 / 0.195

RT

## reward space RV-OV
# best only
S2BWRTmod1bBEST <- lmer(logRT~ cgoalvAvRemBid+cAvBid+(cgoalvAvRemBid+cAvBid|subj_idx), data = a1fMRI[a1fMRI$isChooseBestTrial==1,],
                                 REML=FALSE)

sv1_max <- svd(getME(S2BWRTmod1bBEST, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

#worst only
S2BWRTmod1bWORST <- lmer(logRT~ cgoalvAvRemBid+cAvBid+(cgoalvAvRemBid+cAvBid|subj_idx), data = a1fMRI[a1fMRI$isChooseBestTrial==0,],
                                 REML=FALSE)

sv1_max <- svd(getME(S2BWRTmod1bWORST, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  
#p.val= "kr", show.df = T,
tab_model(S2BWRTmod1bBEST,S2BWRTmod1bWORST,  transform= NULL ,show.stat = TRUE, string.stat="t",   dv.labels = c("log RT Best", "log RT worst"),pred.labels= c("(Intercept)", "Value Difference", "Overall Value"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT Best				log RT worst
Predictors	Estimates	CI	t	p	Estimates	CI	t	p
(Intercept)	1.51	1.43 – 1.59	35.66	<0.001	1.53	1.45 – 1.61	36.82	<0.001
Value Difference	-0.44	-0.53 – -0.36	-10.10	<0.001	-0.40	-0.50 – -0.29	-7.41	<0.001
Overall Value	-0.39	-0.49 – -0.30	-8.07	<0.001	0.38	0.27 – 0.49	6.83	<0.001
Random Effects
σ²	0.14				0.15
τ₀₀	0.05 _{subj_idx}				0.05 _{subj_idx}
τ₁₁	0.02 _{subj_idx.cgoalvAvRemBid}				0.04 _{subj_idx.cgoalvAvRemBid}
	0.03 _{subj_idx.cAvBid}				0.05 _{subj_idx.cAvBid}
ρ₀₁	-0.44				-0.07
	0.48				-0.47
ICC	0.28				0.26
N	30 _{subj_idx}				30 _{subj_idx}
Observations	2133				2137
Marginal R² / Conditional R²	0.090 / 0.344				0.069 / 0.315

standard analysis without interaction with choice goal

S2BWRTmod0 <- lmer(logRT~ cgoalvAvRemBid+cAvBid+isChooseBestTrial+(cgoalvAvRemBid+isChooseBestTrial|subj_idx), data = a1fMRI,
                                 REML=FALSE)

sv1_max <- svd(getME(S2BWRTmod0, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(S2BWRTmod0,p.val= "kr", show.df = T, transform= NULL ,show.stat = TRUE, string.stat="t",   dv.labels = c("log RT"),pred.labels= c("(Intercept)", "Value Difference", "Overall Value", "Best - Worst Condition"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.52	1.44 – 1.59	38.86	31.00	<0.001
Value Difference	-0.41	-0.50 – -0.32	-8.63	31.00	<0.001
Overall Value	-0.01	-0.06 – 0.05	-0.19	4221.00	0.847
Best - Worst Condition	-0.01	-0.09 – 0.06	-0.31	31.00	0.758
Random Effects
σ²	0.16
τ₀₀ _{subj_idx}	0.04
τ₁₁ _{subj_idx.cgoalvAvRemBid}	0.04
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.04
ρ₀₁	-0.25
	-0.24
ICC	0.26
N _{subj_idx}	30
Observations	4270
Marginal R² / Conditional R²	0.041 / 0.293

adding interaction of choice goal and overall value

S2BWRTmod1b <- lmer(logRT~ cgoalvAvRemBid+cAvBid*isChooseBestTrial+(cgoalvAvRemBid+isChooseBestTrial|subj_idx), data = a1fMRI,
                                 REML=FALSE)

sv1_max <- svd(getME(S2BWRTmod1b, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(S2BWRTmod1b,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("log RT"), pred.labels= c("(Intercept)", "Value Difference", "Overall Value", "Best - Worst Condition", "Overall Value: Best - Worst"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.52	1.44 – 1.59	38.93	31.00	<0.001
Value Difference	-0.43	-0.52 – -0.34	-9.18	31.00	<0.001
Overall Value	-0.01	-0.06 – 0.05	-0.24	4222.00	0.809
Best - Worst Condition	-0.01	-0.09 – 0.06	-0.33	31.00	0.740
Overall Value: Best - Worst	-0.76	-0.86 – -0.65	-14.39	3914.00	<0.001
Random Effects
σ²	0.15
τ₀₀ _{subj_idx}	0.04
τ₁₁ _{subj_idx.cgoalvAvRemBid}	0.04
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.04
ρ₀₁	-0.27
	-0.09
ICC	0.27
N _{subj_idx}	30
Observations	4270
Marginal R² / Conditional R²	0.079 / 0.330

nested model to get by condition estimates while controlling for main effect

S2BWRTmod1bn <- lmer(logRT~ cgoalvAvRemBid+isChooseBestTrial/(cAvBid)+(cgoalvAvRemBid+isChooseBestTrial|subj_idx), data = a1fMRI,
                                  REML=FALSE)

sv1_max <- svd(getME(S2BWRTmod1bn, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1) 

tab_model(S2BWRTmod1bn,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("log RT"), pred.labels= c("(Intercept)", "Best - Worst Condition", "Value Difference", "Worst: Overall Value", "Best: Overall Value"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.52	1.44 – 1.59	38.93	31.00	<0.001
Best - Worst Condition	-0.43	-0.52 – -0.34	-9.18	31.00	<0.001
Value Difference	-0.01	-0.09 – 0.06	-0.33	31.00	0.740
Worst: Overall Value	0.37	0.30 – 0.44	9.97	4103.00	<0.001
Best: Overall Value	-0.38	-0.46 – -0.31	-10.28	4155.00	<0.001
Random Effects
σ²	0.15
τ₀₀ _{subj_idx}	0.04
τ₁₁ _{subj_idx.cgoalvAvRemBid}	0.04
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.04
ρ₀₁	-0.27
	-0.09
ICC	0.27
N _{subj_idx}	30
Observations	4270
Marginal R² / Conditional R²	0.079 / 0.330

Figure 1 C overlay part 2

eff_df <- Effect(c("cAvBid", "isChooseBestTrial"), S2BWRTmod1b)
IA <- as.data.frame(eff_df)
IA$cAvBid <- (IA$cAvBid - min(IA$cAvBid))*10

IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
pASVBW <- ggplot(data=IA, aes(x=cAvBid, y=fit , color=isChooseBestTrial )) + geom_line()+scale_colour_manual(name="Goal",values=c("#CC0000","#000000")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cAvBid, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous(limits=c(0.8, 2.1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Overal value") + ylab(" ") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+ theme(legend.justification=c(1,1), legend.position=c(0.95,0.25) )

eff_df <- Effect(c("cgoalvAvRemBid", "isChooseBestTrial"), S2BWRTmod1b, xlevels=list(cgoalvAvRemBid=seq(min(a1fMRI$cgoalvAvRemBid), max(a1fMRI$cgoalvAvRemBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$cgoalvAvRemBid <- (IA$cgoalvAvRemBid - min(IA$cgoalvAvRemBid))*10
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))


## with CI instead of se
pVDBW <- ggplot(data=IA, aes(x=cgoalvAvRemBid, y=fit , color=isChooseBestTrial )) + geom_line()+scale_colour_manual(name="Goal",values=c("#CC0000","#000000")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cgoalvAvRemBid, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous(limits=c(0.8, 2.1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Relative Value") + ylab("log RT") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+ theme(legend.position="none" )

multiplot( pVDBW,pASVBW, cols=2)

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/BASBWfMRI/best_worst_within_allRS.pdf", width = 8, height = 4)#, units = 'cm', res = 200, compression = 'lzw'
multiplot( pVDBW,pASVBW, cols=2)
dev.off()

pdf 
  2

goal value instead of reward value

S2BWRTmod1g <- lmer(logRT~ cgoalvAvRemBid+cavGoalValue+isChooseBestTrial+(cgoalvAvRemBid+isChooseBestTrial|subj_idx), data = a1fMRI,
                                 REML=FALSE)

sv1_max <- svd(getME(S2BWRTmod1g, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(S2BWRTmod1g,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("log RT"), pred.labels= c("(Intercept)", "Value Difference", "Overall Goal Value", "Best - Worst Condition"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.52	1.44 – 1.59	38.98	31.00	<0.001
Value Difference	-0.43	-0.52 – -0.34	-9.19	31.00	<0.001
Overall Goal Value	-0.38	-0.43 – -0.33	-14.39	3912.00	<0.001
Best - Worst Condition	-0.03	-0.11 – 0.04	-0.85	31.00	0.402
Random Effects
σ²	0.15
τ₀₀ _{subj_idx}	0.04
τ₁₁ _{subj_idx.cgoalvAvRemBid}	0.04
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.04
ρ₀₁	-0.27
	-0.09
ICC	0.27
N _{subj_idx}	30
Observations	4270
Marginal R² / Conditional R²	0.079 / 0.329

eff_df <- Effect(c("cavGoalValue", "isChooseBestTrial"), S2BWRTmod1g)
IA <- as.data.frame(eff_df)
IA$cavGoalValue <-  (IA$cavGoalValue - min(IA$cavGoalValue))*10
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
pASVBW <- ggplot(data=IA, aes(x=cavGoalValue, y=fit , color=isChooseBestTrial )) + geom_line()+scale_colour_manual(name="Goal",values=c("#CC0000","#000000")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cavGoalValue, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous(limits=c(0.8, 2.1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Overall value") + ylab(" ") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.justification=c(1,1), legend.position=c(0.95,0.25) )

eff_df <- Effect(c("cgoalvAvRemBid", "isChooseBestTrial"), S2BWRTmod1g, xlevels=list(cgoalvAvRemBid=seq(min(a1fMRI$cgoalvAvRemBid), max(a1fMRI$cgoalvAvRemBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$cgoalvAvRemBid <- (IA$cgoalvAvRemBid - min(IA$cgoalvAvRemBid))*10
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
pVDBW <- ggplot(data=IA, aes(x=cgoalvAvRemBid, y=fit , color=isChooseBestTrial )) + geom_line()+scale_colour_manual(name="Goal",values=c("#CC0000","#000000")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cgoalvAvRemBid, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous(limits=c(0.8, 2.1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Relative Value") + ylab("log RT") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="none")

multiplot(pVDBW, pASVBW, cols=2)

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/BASBWfMRI/best_worst_within_allGS.pdf", width = 8, height = 4)#, units = 'cm', res = 200, compression = 'lzw'
multiplot(pVDBW, pASVBW, cols=2)
dev.off()

pdf 
  2

testing for residual OV effects

S2BWRTmod1allin <- lmer(logRT~ cgoalvAvRemBid+cavGoalValue+ cAvBid+isChooseBestTrial+((cgoalvAvRemBid)+isChooseBestTrial|subj_idx), data = a1fMRI,
                                     REML=FALSE)

sv1_max <- svd(getME(S2BWRTmod1allin, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(S2BWRTmod1allin,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("log RT"), pred.labels= c("(Intercept)", "Value Difference", "Overall Goal Value", "Overall Reward Value", "Best - Worst Condition"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.52	1.44 – 1.59	38.93	31.00	<0.001
Value Difference	-0.43	-0.52 – -0.34	-9.18	31.00	<0.001
Overall Goal Value	-0.38	-0.43 – -0.33	-14.39	3914.00	<0.001
Overall Reward Value	-0.01	-0.06 – 0.05	-0.24	4222.00	0.809
Best - Worst Condition	-0.03	-0.11 – 0.04	-0.85	31.00	0.402
Random Effects
σ²	0.15
τ₀₀ _{subj_idx}	0.04
τ₁₁ _{subj_idx.cgoalvAvRemBid}	0.04
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.04
ρ₀₁	-0.27
	-0.09
ICC	0.27
N _{subj_idx}	30
Observations	4270
Marginal R² / Conditional R²	0.079 / 0.330

model comparison

Table 2 part 2

S2BWRTmodbase <- lmer(logRT~ cgoalvAvRemBid+isChooseBestTrial+(cgoalvAvRemBid+isChooseBestTrial|subj_idx), data = a1fMRI,
                                     REML=FALSE)

sv1_max <- svd(getME(S2BWRTmodbase, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(S2BWRTmodbase,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("log RT"), pred.labels= c("(Intercept)", "Value Difference", "Best - Worst Condition"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.52	1.44 – 1.59	38.89	31.00	<0.001
Value Difference	-0.41	-0.50 – -0.32	-8.65	31.00	<0.001
Best - Worst Condition	-0.01	-0.09 – 0.06	-0.31	31.00	0.757
Random Effects
σ²	0.16
τ₀₀ _{subj_idx}	0.04
τ₁₁ _{subj_idx.cgoalvAvRemBid}	0.04
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.04
ρ₀₁	-0.25
	-0.24
ICC	0.26
N _{subj_idx}	30
Observations	4270
Marginal R² / Conditional R²	0.041 / 0.292

# model comparisons
xxx <- anova(S2BWRTmod0, S2BWRTmod1b, S2BWRTmod1g, S2BWRTmodbase)
row.names(xxx) <- c("VD + C<sub>b-w</sub> (baseline)", "baseline + OV<sub>reward</sub>", "baseline + OV<sub>goal</sub>", "baseline + C<sub>b-w</sub>: OV<sub>reward</sub>")
attributes(xxx)$heading <- NULL
xxx2 <- xxx#[, c(1,2, 6, 7, 8)]
xxx2$AIC <- round(xxx2$AIC)
xxx2$BIC <- round(xxx2$BIC)
xxx2$Chisq <- round(xxx2$Chisq,2)
xxx2[,8] <- round(xxx2[,8],3)
xxx2$logLik <- xxx2$AIC - c(NaN, xxx2$AIC[1:length(xxx2$AIC)-1]) 
allR2s <- rsquared(list(S2BWRTmodbase,S2BWRTmod0, S2BWRTmod1b, S2BWRTmod1g))
xxx2[,1] <- round(allR2s$Conditional,2)
xxx2 <- xxx2[, c(1:4, 6, 8)]
colnames(xxx2) <- c("R<sup>2</sub>", "AIC", "BIC", "dAIC", "Chi<sup>2</sub>", "p")
htmlTable(xxx2)

	R²	AIC	BIC	dAIC	Chi²	p
VD + C_b-w (baseline)	0.19	4396	4460
baseline + OV_reward	0.19	4398	4468	2	0.04	0.847
baseline + OV_goal	0.25	4195	4265	-203	203.31	0
baseline + C_b-w: OV_reward	0.25	4197	4273	2	0.06	0.809

Table 1:

# top
tab_model(BWRTmod1b,S2BWRTmod1b,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("Study 1", "Study 2"), pred.labels= c("(Intercept)", "Value Difference", "Overall Value", "Best - Worst Condition", "Overall Value: Best - Worst"), col.order = c("est", "se","ci",  "stat","df", "p"))

	Study 1					Study 2
Predictors	Estimates	CI	t	df	p	Estimates	CI	t	df	p
(Intercept)	1.42	1.36 – 1.49	41.93	31.00	<0.001	1.52	1.44 – 1.59	38.93	31.00	<0.001
Value Difference	-0.49	-0.55 – -0.43	-16.48	3430.00	<0.001	-0.43	-0.52 – -0.34	-9.18	31.00	<0.001
Overall Value	0.05	-0.01 – 0.11	1.61	3448.00	0.108	-0.01	-0.06 – 0.05	-0.24	4222.00	0.809
Best - Worst Condition	-0.07	-0.14 – 0.00	-1.90	31.00	0.066	-0.01	-0.09 – 0.06	-0.33	31.00	0.740
Overall Value: Best - Worst	-0.75	-0.86 – -0.63	-12.66	3296.00	<0.001	-0.76	-0.86 – -0.65	-14.39	3914.00	<0.001
Random Effects
σ²	0.17					0.15
τ₀₀	0.03 _{subj_idx}					0.04 _{subj_idx}
τ₁₁	0.03 _{subj_idx.isChooseBestTrial2-1}					0.04 _{subj_idx.cgoalvAvRemBid}
						0.04 _{subj_idx.isChooseBestTrial2-1}
ρ₀₁	-0.00 _{subj_idx}					-0.27
						-0.09
ICC	0.18					0.27
N	30 _{subj_idx}					30 _{subj_idx}
Observations	3470					4270
Marginal R² / Conditional R²	0.096 / 0.262					0.079 / 0.330

# bottom
tab_model(BWRTmod1allin, S2BWRTmod1allin,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("Study 1", "Study 2"), pred.labels= c("(Intercept)", "Value Difference", "Overall Goal Value", "Overall Reward Value", "Best - Worst Condition"), col.order = c("est", "se","ci",  "stat","df", "p"))

	Study 1					Study 2
Predictors	Estimates	CI	t	df	p	Estimates	CI	t	df	p
(Intercept)	1.42	1.36 – 1.49	41.97	31.00	<0.001	1.52	1.44 – 1.59	38.93	31.00	<0.001
Value Difference	-0.49	-0.55 – -0.43	-16.48	3430.00	<0.001	-0.43	-0.52 – -0.34	-9.18	31.00	<0.001
Overall Goal Value	-0.37	-0.43 – -0.32	-12.66	3296.00	<0.001	-0.38	-0.43 – -0.33	-14.39	3914.00	<0.001
Overall Reward Value	0.05	-0.01 – 0.11	1.61	3448.00	0.108	-0.01	-0.06 – 0.05	-0.24	4222.00	0.809
Best - Worst Condition	-0.08	-0.15 – -0.01	-2.21	31.00	0.035	-0.03	-0.11 – 0.04	-0.85	31.00	0.402
Random Effects
σ²	0.17					0.15
τ₀₀	0.03 _{subj_idx}					0.04 _{subj_idx}
τ₁₁	0.03 _{subj_idx.isChooseBestTrial2-1}					0.04 _{subj_idx.cgoalvAvRemBid}
						0.04 _{subj_idx.isChooseBestTrial2-1}
ρ₀₁	-0.00 _{subj_idx}					-0.27
						-0.09
ICC	0.18					0.27
N	30 _{subj_idx}					30 _{subj_idx}
Observations	3470					4270
Marginal R² / Conditional R²	0.096 / 0.262					0.079 / 0.330

Liking

As in Study 2, Liking is explained by reward value, not goal value

a1fMRI$LikingOrdinal <- as.factor(a1fMRI$Liking)
#ordered(a1fMRI$LikingOrdinal, levels=c("1", "2", "3", "4", "5"))
a1fMRI$cMaxV <- scale(a1fMRI$MaxBid/10, scale=FALSE, center=TRUE)
fmm1Lik <- clmm(LikingOrdinal ~ isChooseBestTrial+crelGoalvAvRem +cAvBid +     cavGoalValue +cRT+(isChooseBestTrial+cRT|subj_idx), data= a1fMRI[!a1fMRI$Choice1==-1,])

tab_model(fmm1Lik, transform= NULL,show.stat = TRUE, string.stat = "z", pred.labels= c("(Intercept: 1|2)","(Intercept: 2|3)", "(Intercept: 3|4)", "(Intercept: 4|5)","Best - Worst Condition", "Value Difference","Overall Reward Value", "Overall Goal Value",  "RT"))

	LikingOrdinal
Predictors	Log-Odds	CI	z	p
(Intercept: 1\|2)	-2.35	-2.76 – -1.95	-11.33	<0.001
(Intercept: 2\|3)	0.02	-0.38 – 0.42	0.08	0.935
(Intercept: 3\|4)	1.62	1.22 – 2.02	7.89	<0.001
(Intercept: 4\|5)	3.35	2.94 – 3.77	15.81	<0.001
Best - Worst Condition	0.07	-0.11 – 0.25	0.78	0.435
Value Difference	0.60	0.34 – 0.86	4.58	<0.001
Overall Reward Value	5.68	5.36 – 6.00	34.84	<0.001
Overall Goal Value	0.12	-0.14 – 0.39	0.90	0.367
RT	-0.04	-0.22 – 0.15	-0.39	0.693
N _{subj_idx}	29
Observations	4126

fmm1Likb <- clmm(LikingOrdinal ~ isChooseBestTrial+cMaxV +cAvBid +     cavGoalValue +cRT+(isChooseBestTrial+cRT|subj_idx), data= a1fMRI[!a1fMRI$Choice1==-1,])

tab_model(fmm1Likb, transform= NULL,show.stat = TRUE, string.stat = "z", pred.labels= c("(Intercept: 1|2)","(Intercept: 2|3)", "(Intercept: 3|4)", "(Intercept: 4|5)","Best - Worst Condition", "Max Value","Overall Reward Value", "Overall Goal Value",  "RT"))

	LikingOrdinal
Predictors	Log-Odds	CI	z	p
(Intercept: 1\|2)	-2.36	-2.76 – -1.95	-11.32	<0.001
(Intercept: 2\|3)	0.02	-0.38 – 0.42	0.11	0.914
(Intercept: 3\|4)	1.62	1.22 – 2.02	7.88	<0.001
(Intercept: 4\|5)	3.35	2.93 – 3.77	15.77	<0.001
Best - Worst Condition	0.06	-0.12 – 0.24	0.70	0.485
Max Value	0.86	0.49 – 1.22	4.61	<0.001
Overall Reward Value	4.87	4.41 – 5.33	20.74	<0.001
Overall Goal Value	0.09	-0.17 – 0.36	0.68	0.497
RT	-0.04	-0.22 – 0.15	-0.41	0.678
N _{subj_idx}	29
Observations	4126

fMRI analyses

1) Bartra ROI

make covariance table for predictors:

Some basic checks: Predictors are uncorrelated

Figure S4

library(corrplot)
covMatPred <- cor(cbind(a1fMRI$cAvBid, a1fMRI$ChosenvAvRemV, a1fMRI$cavGoalValue, a1fMRI$crelGoalvAvRem), method = "spearman")

covMatPlot <- corrplot(covMatPred, type = "upper", order = "original", diag = FALSE,addCoef.col=TRUE,
         tl.col = "black", tl.srt = 45)

pOVGOVRHist <- ggplot(data=a1fMRI, aes(x=avGoalValue, y= AvBid))+ geom_point(shape="O", color="#C0C0C0")+theme_bw(12)+ 
  xlab("Overall Goal Value") + ylab("Overall Reward Value") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")

ggMarginal(pOVGOVRHist, type = "histogram", fill="#C0C0C0", xparams = list(size=0.1), yparams = list(size=0.1))

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/BASBWfMRI/PredictorScatterOV.pdf", width = 4, height = 4)#, units = 'cm', res = 200, compression = 'lzw'
ggMarginal(pOVGOVRHist, type = "histogram", fill="#C0C0C0", xparams = list(size=0.1), yparams = list(size=0.1))
dev.off()

pdf 
  2

pRVGRVRHist <- ggplot(data=a1fMRI, aes(x=relGoalvAvRem, y= ChosenvAvRemV))+ geom_point(shape="O", color="#C0C0C0")+theme_bw(12)+ 
  xlab("Chosen vs Unchosen Goal Value") + ylab("Chosen vs Unchosen Reward Value") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")

ggMarginal(pRVGRVRHist, type = "histogram", fill="#C0C0C0", xparams = list(size=0.1), yparams = list(size=0.1))

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/BASBWfMRI/PredictorScatterRV.pdf", width = 4, height = 4)#, units = 'cm', res = 200, compression = 'lzw'
ggMarginal(pRVGRVRHist, type = "histogram", fill="#C0C0C0", xparams = list(size=0.1), yparams = list(size=0.1))
dev.off()

pdf 
  2

Testing whether reward or goal values (or both) explain value network (BOLD) activity

Baseline model

# baseline model
BWBACallInModctestbase <- lmer(binConjunc_PvNxDECxRECxMONxPRI~ isChooseBestTrial+cRT+(isChooseBestTrial+cRT|subj_idx), data = a1fMRI,
                                           REML=FALSE)

sv1_max <- svd(getME(BWBACallInModctestbase, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(BWBACallInModctestbase,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("BOLD Activity Value Network"), pred.labels= c("(Intercept)", "Best - Worst Condition" , "RT"), col.order = c("est", "se","ci",  "stat","df", "p"))

	BOLD Activity Value Network
Predictors	Estimates	CI	t	df	p
(Intercept)	0.00	-0.08 – 0.08	0.02	31.00	0.985
Best - Worst Condition	-0.02	-0.11 – 0.07	-0.52	31.00	0.608
RT	-0.03	-0.13 – 0.07	-0.52	30.00	0.607
Random Effects
σ²	0.94
τ₀₀ _{subj_idx}	0.04
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.03
τ₁₁ _{subj_idx.cRT}	0.03
ρ₀₁	0.17
	-0.08
ICC	0.06
N _{subj_idx}	30
Observations	4270
Marginal R² / Conditional R²	0.000 / 0.059

Reward Value model

# Reward Value model:
BWBACallInModctest0 <- lmer(binConjunc_PvNxDECxRECxMONxPRI~ isChooseBestTrial+cAvBid+ ChosenvAvRemV +cRT+(isChooseBestTrial+cRT|subj_idx), data = a1fMRI,
                                      REML=FALSE)

sv1_max <- svd(getME(BWBACallInModctest0, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(BWBACallInModctest0,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("BOLD Activity Value Network"), pred.labels= c("(Intercept)", "Best - Worst Condition" , "Overall Reward Value", "Relative Reward Value", "RT"), col.order = c("est", "se","ci",  "stat","df", "p"))

	BOLD Activity Value Network
Predictors	Estimates	CI	t	df	p
(Intercept)	0.00	-0.08 – 0.08	0.01	31.00	0.993
Best - Worst Condition	-0.02	-0.11 – 0.08	-0.36	41.00	0.720
Overall Reward Value	0.21	0.09 – 0.34	3.26	3861.00	0.001
Relative Reward Value	-0.02	-0.15 – 0.10	-0.36	4233.00	0.718
RT	-0.03	-0.13 – 0.07	-0.52	30.00	0.605
Random Effects
σ²	0.94
τ₀₀ _{subj_idx}	0.04
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.03
τ₁₁ _{subj_idx.cRT}	0.03
ρ₀₁	0.15
	-0.05
ICC	0.06
N _{subj_idx}	30
Observations	4270
Marginal R² / Conditional R²	0.003 / 0.062

Goal Value Model

# Goal value model
BWBACallInModctest2 <- lmer(binConjunc_PvNxDECxRECxMONxPRI~ isChooseBestTrial+cavGoalValue+ crelGoalvAvRem +cRT+(isChooseBestTrial+cRT|subj_idx), data = a1fMRI,
                                           REML=FALSE)

sv1_max <- svd(getME(BWBACallInModctest2, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(BWBACallInModctest2,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("BOLD Activity Value Network"), pred.labels= c("(Intercept)", "Best - Worst Condition" , "Overall Goal Value", "Relative Goal Value", "RT"), col.order = c("est", "se","ci",  "stat","df", "p"))

	BOLD Activity Value Network
Predictors	Estimates	CI	t	df	p
(Intercept)	0.00	-0.08 – 0.08	0.02	31.00	0.984
Best - Worst Condition	-0.01	-0.10 – 0.07	-0.31	32.00	0.758
Overall Goal Value	0.16	0.03 – 0.29	2.41	1635.00	0.016
Relative Goal Value	0.23	0.10 – 0.36	3.49	4167.00	<0.001
RT	0.02	-0.08 – 0.12	0.39	32.00	0.698
Random Effects
σ²	0.94
τ₀₀ _{subj_idx}	0.04
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.03
τ₁₁ _{subj_idx.cRT}	0.03
ρ₀₁	0.15
	-0.06
ICC	0.06
N _{subj_idx}	30
Observations	4270
Marginal R² / Conditional R²	0.004 / 0.060

Test for residual Overall Reward Value effects

# goal values plus OVr
BWBACallInModctestOVg <- lmer(binConjunc_PvNxDECxRECxMONxPRI~ isChooseBestTrial+cAvBid+cavGoalValue+ crelGoalvAvRem+cRT+(isChooseBestTrial+cRT|subj_idx), data = a1fMRI,REML=FALSE)

sv1_max <- svd(getME(BWBACallInModctestOVg, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(BWBACallInModctestOVg,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("BOLD Activity Value Network"), pred.labels= c("(Intercept)", "Best - Worst Condition" ,"Overall Reward Value", "Overall Goal Value", "Relative Goal Value", "RT"), col.order = c("est", "se","ci",  "stat","df", "p"))

	BOLD Activity Value Network
Predictors	Estimates	CI	t	df	p
(Intercept)	0.00	-0.08 – 0.08	0.01	31.00	0.991
Best - Worst Condition	-0.01	-0.10 – 0.07	-0.31	32.00	0.755
Overall Reward Value	0.21	0.08 – 0.34	3.18	3845.00	0.001
Overall Goal Value	0.16	0.03 – 0.29	2.43	1625.00	0.015
Relative Goal Value	0.22	0.09 – 0.36	3.37	4173.00	0.001
RT	0.02	-0.08 – 0.12	0.37	32.00	0.713
Random Effects
σ²	0.93
τ₀₀ _{subj_idx}	0.04
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.03
τ₁₁ _{subj_idx.cRT}	0.03
ρ₀₁	0.14
	-0.03
ICC	0.06
N _{subj_idx}	30
Observations	4270
Marginal R² / Conditional R²	0.007 / 0.063

Full model with all predictors

Table 4

# full model w all predictors
BWBACallInModc <- lmer(binConjunc_PvNxDECxRECxMONxPRI~ isChooseBestTrial+cAvBid+ ChosenvAvRemV +cavGoalValue+ crelGoalvAvRem+cRT+(isChooseBestTrial+cRT|subj_idx), data = a1fMRI,
                                      REML=FALSE)

sv1_max <- svd(getME(BWBACallInModc, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(BWBACallInModc,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("BOLD Activity Value Network"), pred.labels= c("(Intercept)", "Best - Worst Condition" ,"Overall Reward Value", "Relative Reward Value", "Overall Goal Value", "Relative Goal Value", "RT"), col.order = c("est", "se","ci",  "stat","df", "p"))

	BOLD Activity Value Network
Predictors	Estimates	CI	t	df	p
(Intercept)	0.00	-0.08 – 0.08	0.01	31.00	0.992
Best - Worst Condition	-0.01	-0.10 – 0.09	-0.13	42.00	0.896
Overall Reward Value	0.21	0.08 – 0.34	3.15	3844.00	0.002
Relative Reward Value	-0.03	-0.16 – 0.10	-0.44	4225.00	0.657
Overall Goal Value	0.16	0.03 – 0.29	2.44	1632.00	0.015
Relative Goal Value	0.22	0.09 – 0.36	3.37	4173.00	0.001
RT	0.02	-0.08 – 0.12	0.37	32.00	0.715
Random Effects
σ²	0.93
τ₀₀ _{subj_idx}	0.04
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.03
τ₁₁ _{subj_idx.cRT}	0.03
ρ₀₁	0.15
	-0.03
ICC	0.06
N _{subj_idx}	30
Observations	4270
Marginal R² / Conditional R²	0.007 / 0.063

Model Comparison

The model with Relative and Overall Goal Value and Overall Reward Value performs best.

Table 3

xxx <- anova(BWBACallInModctest0, BWBACallInModctestbase,BWBACallInModctestOVg,BWBACallInModctest2, BWBACallInModc)
row.names(xxx) <- c("C<sub>b-w</sub> + RT (baseline)", "baseline + RV<sub>reward</sub> + OV<sub>reward</sub>", "baseline+ RV<sub>goal</sub> + OV<sub>goal</sub>", "baseline+ RV<sub>goal</sub> + OV<sub>goal</sub> + OV<sub>reward</sub>", "baseline+ RV<sub>goal</sub> + OV<sub>goal</sub> +RV<sub>reward</sub> + OV<sub>reward</sub>")
attributes(xxx)$heading <- NULL
xxx2 <- xxx#[, c(1,2, 6, 7, 8)]
xxx2$AIC <- round(xxx2$AIC)
xxx2$BIC <- round(xxx2$BIC)
xxx2$Chisq <- round(xxx2$Chisq,2)
xxx2[,8] <- round(xxx2[,8],3)
xxx2$logLik <- xxx2$AIC - c(NaN, xxx2$AIC[1:length(xxx2$AIC)-1]) 
allR2s <- rsquared(list(BWBACallInModctestbase,BWBACallInModctest0,BWBACallInModctest2, BWBACallInModctestOVg, BWBACallInModc))
xxx2[,1] <- round(allR2s$Conditional,3)
xxx2 <- xxx2[, c(1:4, 6, 8)]
colnames(xxx2) <- c("R<sup>2</sup>", "AIC", "BIC", "dAIC", "Chi<sup>2</sup>", "p")
htmlTable(xxx2)

	R²	AIC	BIC	dAIC	Chi²	p
C_b-w + RT (baseline)	0.062	11976	12040
baseline + RV_reward + OV_reward	0.066	11969	12045	-7	11	0.004
baseline+ RV_goal + OV_goal	0.065	11963	12040	-6	5.68	0
baseline+ RV_goal + OV_goal + OV_reward	0.069	11955	12038	-8	10.2	0.001
baseline+ RV_goal + OV_goal +RV_reward + OV_reward	0.069	11957	12046	2	0.2	0.657

Figure 2

# refit relevant model with lmerTest attached to get the p-values... 

library(lmerTest)
BWBACallInModc <- lmer(binConjunc_PvNxDECxRECxMONxPRI~ isChooseBestTrial+cAvBid+ ChosenvAvRemV +cavGoalValue+ crelGoalvAvRem+cRT+(isChooseBestTrial+cRT|subj_idx), data = a1fMRI,
                                      REML=FALSE)

TEST <- as.data.frame(coef(summary(BWBACallInModc)))
TEST$effects <- rownames(TEST)
TEST$se <- TEST[,2]
TEST$effects<- revalue(TEST$effects, c("(Intercept)"="Intercept", "isChooseBestTrial2-1"= "Task goal",  "cAvBid"="OV[reward]", "ChosenvAvRemV"="RV[reward]",  "cavGoalValue"= "OV[goal]", "crelGoalvAvRem"= "RV[goal]"))
TEST$effects<- ordered(TEST$effects, levels=c("Intercept", "Task goal","RV[reward]",   "OV[reward]", "RV[goal]", "OV[goal]", "cRT"))
TEST <- TEST[c(3:6), ]
TEST$effects <- factor(TEST$effects) 
TEST$isSig <- rep("*", length(TEST$Estimate))
TEST$isSig[TEST[,5]>0.05] <- NA_real_ 
TEST$isSig2 <- TEST$isSig
TEST$isSig2[TEST[,5]<0.01]<- "**" 

EFFPLOT<- ggplot(data=TEST, aes(x=effects, y=Estimate, fill=effects, group=effects, color=effects, shape=effects)) + geom_point(size=5,position=position_dodge(0.6))+theme_bw(12)+ scale_color_manual(values=c("#70AD47", "#70AD47","#313A90",  "#313A90"))+
geom_errorbar(aes(max = Estimate + se, min = Estimate -se), position=position_dodge(0.6),width=0.1)+scale_shape_manual(values = c(19,15,  19, 15))+  
  xlab("") + ylab("Regression Coefficient") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position = "none")+
  geom_hline(aes(yintercept=0), size=0.2,linetype=2, color="black")+
  geom_text(aes(label=isSig2, y=0.35),position=position_dodge(0.6), size=6,show.legend = F) 

EFFPLOT

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/BASBWfMRI/BARTRAROIEFFECTS_chosen_vs_unchosen.pdf", width = 5, height = 3)
EFFPLOT
dev.off()

pdf 
  2

detach("package:lmerTest", unload = TRUE) # detaching it again for tables to work.

2) Follow up analyses

vmPFC and ventral Striatum separately

BWvMPFCCallInMod <- lmer(BARTRA_vmPFC~ isChooseBestTrial+cAvBid+ ChosenvAvRemV +cavGoalValue+ crelGoalvAvRem+cRT+(isChooseBestTrial+cRT|subj_idx), data = a1fMRI,
                                        REML=FALSE)

sv1_max <- svd(getME(BWvMPFCCallInMod, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

library(lmerTest)
BWvMPFCCallInModplot <- lmer(BARTRA_vmPFC~ isChooseBestTrial+cAvBid+ ChosenvAvRemV +cavGoalValue+ crelGoalvAvRem+cRT+(isChooseBestTrial+cRT|subj_idx), data = a1fMRI,
                                        REML=FALSE)

TEST <- as.data.frame(coef(summary(BWvMPFCCallInModplot)))
TEST$effects <- rownames(TEST)
TEST$se <- TEST[,2]
TEST$effects<- revalue(TEST$effects, c("(Intercept)"="Intercept", "isChooseBestTrial2-1"= "Task goal",  "cAvBid"="OVr", "ChosenvAvRemV"="RVr",  "cavGoalValue"= "OVg", "crelGoalvAvRem"= "RVg"))
TEST$effects<- ordered(TEST$effects, levels=c("Intercept", "Task goal","RVr",  "RVg", "OVr",  "OVg", "cRT"))
TEST <- TEST[c(3:6), ]
TEST$effects <- factor(TEST$effects) 

TESTvmPFC <- TEST

detach("package:lmerTest", unload = TRUE)

Table S3

BWStrCallInMod <- lmer(BARTRA_NAcc~ isChooseBestTrial+cAvBid+ ChosenvAvRemV +cavGoalValue+ crelGoalvAvRem+cRT+(isChooseBestTrial+cRT|subj_idx), data = a1fMRI,
                                      REML=FALSE)

sv1_max <- svd(getME(BWStrCallInMod, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(BWvMPFCCallInMod, BWStrCallInMod,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("BOLD Activity vmPFC", "BOLD Activity vStr"), pred.labels= c("(Intercept)", "Best - Worst Condition" ,"Overall Reward Value", "Relative Reward Value", "Overall Goal Value", "Relative Goal Value", "RT"), col.order = c("est", "se","ci",  "stat","df", "p"))

	BOLD Activity vmPFC					BOLD Activity vStr
Predictors	Estimates	CI	t	df	p	Estimates	CI	t	df	p
(Intercept)	-0.56	-0.70 – -0.42	-7.76	31.00	<0.001	-0.41	-0.55 – -0.27	-5.73	31.00	<0.001
Best - Worst Condition	-0.06	-0.21 – 0.09	-0.78	46.00	0.439	0.08	-0.10 – 0.27	0.90	37.00	0.374
Overall Reward Value	0.31	0.07 – 0.55	2.52	3747.00	0.012	0.28	0.08 – 0.48	2.80	4022.00	0.005
Relative Reward Value	0.05	-0.18 – 0.28	0.43	4202.00	0.669	-0.13	-0.32 – 0.07	-1.28	4237.00	0.200
Overall Goal Value	0.33	0.09 – 0.56	2.74	972.00	0.006	0.08	-0.12 – 0.28	0.76	2633.00	0.445
Relative Goal Value	0.33	0.09 – 0.57	2.70	4195.00	0.007	0.30	0.10 – 0.50	2.99	4133.00	0.003
RT	0.20	0.01 – 0.39	2.08	32.00	0.045	-0.17	-0.33 – -0.01	-2.12	32.00	0.042
Random Effects
σ²	3.20					2.15
τ₀₀	0.13 _{subj_idx}					0.13 _{subj_idx}
τ₁₁	0.05 _{subj_idx.isChooseBestTrial2-1}					0.17 _{subj_idx.isChooseBestTrial2-1}
	0.11 _{subj_idx.cRT}					0.08 _{subj_idx.cRT}
ρ₀₁	0.28					0.01
	-0.09					0.09
ICC	0.05					0.08
N	30 _{subj_idx}					30 _{subj_idx}
Observations	4270					4270
Marginal R² / Conditional R²	0.007 / 0.056					0.009 / 0.088

library(lmerTest)
BWStrCallInModplot <- lmer(BARTRA_NAcc~ isChooseBestTrial+cAvBid+ ChosenvAvRemV +cavGoalValue+ crelGoalvAvRem+cRT+(isChooseBestTrial+cRT|subj_idx), data = a1fMRI,
                                      REML=FALSE)

TEST <- as.data.frame(coef(summary(BWStrCallInModplot)))
TEST$effects <- rownames(TEST)
TEST$se <- TEST[,2]
TEST$effects<- revalue(TEST$effects, c("(Intercept)"="Intercept", "isChooseBestTrial2-1"= "Task goal",  "cAvBid"="OVr", "ChosenvAvRemV"="RVr",  "cavGoalValue"= "OVg", "crelGoalvAvRem"= "RVg"))
TEST$effects<- ordered(TEST$effects, levels=c("Intercept", "Task goal","RVr",  "RVg", "OVr",  "OVg", "cRT"))
TEST <- TEST[c(3:6), ]
TEST$effects <- factor(TEST$effects) # get rid of levels I don't need
#TEST$effects<- ordered(TEST$effects, levels=c("Intercept", "Task goal","RVr",  "RVg", "OVr",  "OVg", "cRT"))

TEST_vStr <- TEST

Figure S5

TESTvmPFC$Region <- rep("vmPFC", length(TESTvmPFC$Estimate))
TEST_vStr$Region <- rep("vStr", length(TEST_vStr$Estimate))

TEST2 <- rbind(TESTvmPFC, TEST_vStr)
TEST2$effects<- ordered(TEST2$effects, levels=c("Intercept", "Task goal","RVr",   "OVr", "RVg", "OVg", "cRT"))
TEST2$isSig <- rep("*", length(TEST2$Estimate))
#TEST2$isSig[TEST2$effects=="RVg"] <- 0.6
#TEST2$isSig[TEST2$effects=="OVg"] <- 0.5
TEST2$Loc <- rep(0.55, length(TEST2$Estimate))
TEST2$Loc[TEST2$effects=="RVg"] <- 0.6
TEST2$Loc[TEST2$effects=="OVg"] <- 0.5
TEST2$isSig[TEST2[,5]>0.05] <- NA_real_ 
TEST2$isSig2 <- TEST2$isSig
TEST2$isSig2[TEST2[,5]<0.01]<- "**" 
#,position=position_dodge(0.6) position=position_dodge(0.6),geom_point(aes(x=Region, y=isSig), shape=8)+
EFFPLOTIA <- ggplot(data=TEST2, aes(x=Region, y=Estimate, fill=effects, group=effects, color=effects, shape=effects)) + geom_point(size=7,position=position_dodge(0.6))+theme_bw(14)+ scale_color_manual(values=c("#70AD47","#70AD47","#313A90","#313A90","#4a58db",  "#70AD47", "#88d356"))+#facet_wrap(~ssConfidence, nrow=1)+#geom_ribbon(data=IA, aes(x=block.c, max = upper, min = lower,alpha=0.1, fill=sDirection),show.legend =FALSE)
  geom_errorbar(aes(max = Estimate + se, min = Estimate -se), position=position_dodge(0.6),width=0.2)+scale_shape_manual(values = c( 16, 15, 16,15))+ #scale_fill_brewer() + 
  xlab("") + ylab("Regression Coefficient") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + #theme(legend.position="east outside")+
  geom_hline(aes(yintercept=0), size=0.2,linetype=2, color="black")+ #geom_line(data=TEST2, aes(x=Region, y=0), size=0.2,linetype=2, color="black") +
  #theme(axis.text.x = element_text(angle = 60, hjust = 1))+
  geom_text(aes(label=isSig2, y=0.5),position=position_dodge(0.6), size=10,show.legend = F) 
# +  geom_point(aes(x=Region, y=isSig2), shape=8, position= 0.3 )

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/BASBWfMRI/BARTRAROIEFFECTS_chosen_vs_unchosen_all_Regions_Strvs_vmPFC.pdf", width = 6, height = 4)#, units = 'cm', res = 200, compression = 'lzw'
EFFPLOTIA
dev.off()

pdf 
  2

detach("package:lmerTest", unload = TRUE)
EFFPLOTIA

Interaction model

no significant interaction with region when comparing effects in vmPFC and vStr

a2R <- rbind(a1fMRI, a1fMRI)

a2R$Region <- rep(c("vmPFC", "VStr"), each=length(a1fMRI$subj_idx))
a2R$Region <- as.factor(a2R$Region)
a2R$Region <- ordered(a2R$Region, levels= c("vmPFC", "VStr"))
contrasts(a2R$Region) <- contr.poly(2)#contr.sdif(2)
a2R$BOLD <- a2R$BARTRA_vmPFC
a2R$BOLD[a2R$Region=="VStr"] <- a2R$BARTRA_NAcc[a2R$Region=="VStr"] 
contrasts(a2R$isChooseBestTrial) <- contr.sdif(2)

library(lmerTest)
BARTRASubMod0 <- lmer(BOLD~ (isChooseBestTrial+cAvBid+ ChosenvAvRemV +cavGoalValue+ crelGoalvAvRem)*Region+cRT+(isChooseBestTrial+cRT+Region|subj_idx), data = a2R,
                                     REML=FALSE)

sv1_max <- svd(getME(BARTRASubMod0, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)


xxx <- anova(BARTRASubMod0)
row.names(xxx) <- c("Best-Worst Condition", "OV<sub>reward</sub>", "RV<sub>reward</sub>","OV<sub>goal</sub>","RV<sub>goal</sub>", "Region", "RT","C<sub>b-w</sub>: Region", "OV<sub>reward</sub>:Region", "RV<sub>reward</sub>:Region","OV<sub>goal</sub>:Region","RV<sub>goal</sub>:Region")
attributes(xxx)$heading <- NULL
xxx[, c(1:5)] <- round(xxx[, c(1:5)],2)
xxx[, 6] <- round(xxx[, 6],3)
htmlTable(xxx)

	Sum Sq	Mean Sq	NumDF	DenDF	F value	Pr(>F)
Best-Worst Condition	0.05	0.05	1	36.51	0.02	0.897
OV_reward	38.69	38.69	1	7592.14	14.4	0
RV_reward	0.82	0.82	1	8415.16	0.3	0.581
OV_goal	13.88	13.88	1	3314.5	5.17	0.023
RV_goal	41.06	41.06	1	8276.51	15.29	0
Region	11.4	11.4	1	29.76	4.25	0.048
RT	0.01	0.01	1	32.82	0	0.952
C_b-w: Region	7.61	7.61	1	8426.47	2.83	0.092
OV_reward:Region	0.06	0.06	1	4437.32	0.02	0.878
RV_reward:Region	2.52	2.52	1	8439.63	0.94	0.333
OV_goal:Region	5.93	5.93	1	8430.75	2.21	0.137
RV_goal:Region	1.86	1.86	1	8436.61	0.69	0.406

detach("package:lmerTest", unload = TRUE)

BARTRASubMod0t <- lmer(BOLD~ Region/(isChooseBestTrial+cAvBid+ cavGoalValue+ChosenvAvRemV + crelGoalvAvRem)+cRT+(isChooseBestTrial+cRT+Region|subj_idx), data = a2R,
                                      REML=FALSE)

sv1_max <- svd(getME(BARTRASubMod0t, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

#p.val= "kr", show.df = T ,
tab_model(BARTRASubMod0t,show.stat = TRUE, string.stat="t", dv.labels = c("BOLD"), pred.labels= c("(Intercept)", "Region (linear)", "RT", "vmPFC: Best - Worst Condition","vStr: Best - Worst Condition" ,"vmPFC: Overall Reward Value", "vStr: Overall Reward Value","vmPFC:Overall Goal Value","vStr: Overall Goal Value","vmPFC: Relative Reward Value", "vStr: Relative Reward Value",  "vmPFC: Relative Goal Value",  "vStr: Relative Goal Value"), col.order = c("est", "se","ci",  "stat","df", "p"))

	BOLD
Predictors	Estimates	CI	t	p
(Intercept)	-0.48	-0.60 – -0.37	-8.20	<0.001
Region (linear)	0.11	0.01 – 0.21	2.06	0.039
RT	0.00	-0.14 – 0.14	0.06	0.953
vmPFC: Best - Worst Condition	-0.06	-0.22 – 0.10	-0.75	0.455
vStr: Best - Worst Condition	0.08	-0.08 – 0.24	0.97	0.331
vmPFC: Overall Reward Value	0.31	0.09 – 0.53	2.82	0.005
vStr: Overall Reward Value	0.29	0.07 – 0.50	2.61	0.009
vmPFC:Overall Goal Value	0.28	0.08 – 0.49	2.68	0.007
vStr: Overall Goal Value	0.07	-0.13 – 0.28	0.70	0.482
vmPFC: Relative Reward Value	0.03	-0.18 – 0.24	0.29	0.773
vStr: Relative Reward Value	-0.12	-0.33 – 0.10	-1.07	0.284
vmPFC: Relative Goal Value	0.25	0.03 – 0.47	2.24	0.025
vStr: Relative Goal Value	0.38	0.16 – 0.59	3.39	0.001
Random Effects
σ²	2.69
τ₀₀ _{subj_idx}	0.09
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.09
τ₁₁ _{subj_idx.cRT}	0.09
τ₁₁ _{subj_idx.Region.L}	0.07
ρ₀₁	0.14
	-0.02
	-0.09
ICC	0.05
N _{subj_idx}	30
Observations	8540
Marginal R² / Conditional R²	0.007 / 0.058

Gradient analysis Striatum

Table S5

### Data Prep ###
a3 <- rbind(a1fMRI, a1fMRI, a1fMRI)

a3$Region <- rep(c("DC", "VSs", "VSi"), each=length(a1fMRI$subj_idx))
a3$Region <- as.factor(a3$Region)
a3$Region <- ordered(a3$Region, levels= c( "DC", "VSs", "VSi"))
contrasts(a3$Region) <- contr.poly(3)
a3$BOLD <- a3$FCN_Str
a3$BOLD[a3$Region=="VSs"] <- a3$DMN_Str[a3$Region=="VSs"] 
a3$BOLD[a3$Region=="VSi"] <- a3$LN_Str[a3$Region=="VSi"] 
contrasts(a3$isChooseBestTrial) <- contr.sdif(2)

a9 <- rbind(a3, a3)

a9$CATValue <- rep(c("OVr","RVg"), each=length(a3$subj_idx))
a9$CATValue <- as.factor(a9$CATValue)
a9$CATValue <- ordered(a9$CATValue, levels= c( "OVr", "RVg"))
contrasts(a9$CATValue ) <- contr.sdif(2)
a9$ContValue  <- a9$cAvBid 
a9$ContValue[ a9$CATValue =="RVg"] <- a9$crelGoalvAvRem[ a9$CATValue =="RVg"]
contrasts(a9$isChooseBestTrial) <- contr.sdif(2)
contrasts(a9$Region) <- contr.poly(3) # to test for linear and quadratic effects

###################### THIS IS THE GRADIENT ANALYSIS!!!! ################################

library(lmerTest)
Gradient3by3Mod0 <- lmer(BOLD~ Region*CATValue:ContValue+cRT+(cRT|subj_idx), data = a9,
                                        REML=FALSE)

yyy <- summary(Gradient3by3Mod0)
yyy2 <- yyy$coefficients

row.names(yyy2) <-  c("(Intercept)", "Region (linear)","Region (quadratic)","RT", "Overall Reward Value","Relative Goal Value" ,"Region (linear): Overall Reward Value", "Region (quadratic): Overall Reward Value","Region (linear): Relative Goal Value",  "Region (quadratic): Relative Goal Value")
yyy2 <- yyy2[,c(1,2, 4, 3,5)]
colnames(yyy2) <- c("Estimates", "SE","t", "df",  "p")
yyy2[,c(1:4)] <- round(yyy2[,c(1:4)],2)
yyy2[,c(5)] <- round(yyy2[,c(5)],3)
htmlTable(yyy2)

	Estimates	SE	t	df	p
(Intercept)	-0.4	0.06	-6.26	29.61	0
Region (linear)	-0.11	0.02	-6.28	25559.83	0
Region (quadratic)	-0.08	0.02	-4.5	25559.83	0
RT	-0.21	0.07	-3	30.15	0.005
Overall Reward Value	0.35	0.06	5.54	25479.49	0
Relative Goal Value	0.2	0.07	2.99	25531.34	0.003
Region (linear): Overall Reward Value	-0.14	0.1	-1.41	25559.83	0.158
Region (quadratic): Overall Reward Value	-0.21	0.1	-2.1	25559.83	0.036
Region (linear): Relative Goal Value	0.24	0.11	2.19	25559.83	0.029
Region (quadratic): Relative Goal Value	0.02	0.11	0.2	25559.83	0.844

xxx <- anova(Gradient3by3Mod0)
row.names(xxx) <- c("Region", "RT", "Value Type","Region: Value Type")
attributes(xxx)$heading <- NULL
xxx[, c(1:5)] <- round(xxx[, c(1:5)],2)
xxx[, 6] <- round(xxx[, 6],3)
htmlTable(xxx)

	Sum Sq	Mean Sq	NumDF	DenDF	F value	Pr(>F)
Region	170.7	85.35	2	25559.83	29.88	0
RT	25.75	25.75	1	30.15	9.01	0.005
Value Type	112.9	56.45	2	25505.45	19.76	0
Region: Value Type	32.08	8.02	4	25559.83	2.81	0.024

detach("package:lmerTest", unload = TRUE)

Estimates by Regions separately

Dorsal Caudate

library(lmerTest)
BWFCN_StrCallInMod <- lmer(FCN_Str~ isChooseBestTrial+cAvBid+ crelGoalvAvRem+cRT+(isChooseBestTrial|subj_idx), data = a1fMRI,
                                          REML=FALSE)

sv1_max <- svd(getME(BWFCN_StrCallInMod, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

TEST <- as.data.frame(coef(summary(BWFCN_StrCallInMod)))
TEST$effects <- rownames(TEST)
TEST$se <- TEST[,2]
TEST$effects<- revalue(TEST$effects, c("(Intercept)"="Intercept", "isChooseBestTrial2-1"= "Task goal",  "cAvBid"="OVr", "crelGoalvAvRem"= "RVg"))
TEST$effects<-  ordered(TEST$effects, levels=c("Intercept", "Task goal", "RVg", "OVr",  "cRT"))
TEST <- TEST[c(3:4), ]
TEST$effects <- factor(TEST$effects) 
TESTDC <- TEST

Ventral Striatum (superior)

BWDMN_StrCallInMod <- lmer(DMN_Str~ isChooseBestTrial+cAvBid+crelGoalvAvRem+cRT+(isChooseBestTrial|subj_idx), data = a1fMRI,
                                          REML=FALSE)

sv1_max <- svd(getME(BWDMN_StrCallInMod, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)


TEST <- as.data.frame(coef(summary(BWDMN_StrCallInMod)))
TEST$effects <- rownames(TEST)
TEST$se <- TEST[,2]
TEST$effects<- revalue(TEST$effects, c("(Intercept)"="Intercept", "isChooseBestTrial2-1"= "Task goal",  "cAvBid"="OVr", "crelGoalvAvRem"= "RVg"))
TEST$effects<-  ordered(TEST$effects, levels=c("Intercept", "Task goal", "RVg", "OVr",  "cRT"))
TEST <- TEST[c(3:4), ]
TEST$effects <- factor(TEST$effects) 
TESTVSS <- TEST

ventral Striatum (inferior)

BWLN_StrCallInMod <- lmer(LN_Str~ isChooseBestTrial+cAvBid +crelGoalvAvRem+cRT+(isChooseBestTrial|subj_idx), data = a1fMRI, 
                                         REML=FALSE)

sv1_max <- svd(getME(BWLN_StrCallInMod, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)


TEST <- as.data.frame(coef(summary(BWLN_StrCallInMod)))
TEST$effects <- rownames(TEST)
TEST$se <- TEST[,2]
TEST$effects<- revalue(TEST$effects, c("(Intercept)"="Intercept", "isChooseBestTrial2-1"= "Task goal",  "cAvBid"="OVr", "crelGoalvAvRem"= "RVg"))
TEST$effects<-  ordered(TEST$effects, levels=c("Intercept", "Task goal", "RVg", "OVr",  "cRT"))
TEST <- TEST[c(3:4), ]
TEST$effects <- factor(TEST$effects) 


TESTVSI <- TEST

detach("package:lmerTest", unload = TRUE)

FIGURE 3

## integrate all effects into 1 figure
#"DC", "VSs", "VSi"
TESTDC$Region <- rep("Dorsal\n Caudate", length(TESTDC$Estimate))
TESTVSS$Region <- rep("Ventral\n Striatum\n (superior)", length(TESTVSS$Estimate))
TESTVSI$Region <- rep("Ventral\n Striatum\n (inferior)", length(TESTVSI$Estimate))

TEST2 <- rbind(TESTDC, TESTVSS, TESTVSI)
TEST2$Region<- ordered(TEST2$Region, levels=c("Dorsal\n Caudate", "Ventral\n Striatum\n (superior)", "Ventral\n Striatum\n (inferior)"))

TEST3 <- TEST2[!TEST2$effects=="RVr" & !TEST2$effects=="OVg" & !TEST2$effects=="cRT",] 
TEST3$isSig <- rep("*", length(TEST3$Estimate))
TEST3$isSig[TEST3[,5]>0.05] <- NA_real_ 
TEST3$isSig2 <- TEST3$isSig
TEST3$isSig2[TEST3[,5]<0.01]<- "**" 
TEST3$isSig2[TEST3[,5]<0.001]<- "***" 
EFFPLOTIA <- ggplot(data=TEST3, aes(x=Region, y=Estimate, fill=effects, group=effects, color=effects, shape=effects)) + geom_point(size=5)+theme_bw(14)+ scale_color_manual(values=c("#313A90", "#70AD47", "#88d356"))+
  geom_errorbar(aes(max = Estimate + se, min = Estimate -se), width=0.1)+scale_shape_manual(values = c(19, 15, 18))+ 
  xlab("") + ylab("Regression Coefficient") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + 
 geom_hline(yintercept = 0, size=0.2,linetype=2, color="black") +
  geom_text(aes(label=isSig2, y=0.6), size=6,show.legend = F) 

EFFPLOTIA

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/BASBWfMRI/BARTRAROIEFFECTS_chosen_vs_unchosen_all_Regions_Str.pdf", width = 6, height = 4)
EFFPLOTIA
dev.off()

pdf 
  2

3) Set Appraisal

Value network

BWBACallInModLik0 <- lmer(binConjunc_PvNxDECxRECxMONxPRI~ cLiking+(cLiking|subj_idx), data = a1fMRI,
                                         REML=FALSE)

sv1_max <- svd(getME(BWBACallInModLik0, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(BWBACallInModLik0,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("BOLD Activity Value Network"), pred.labels= c("(Intercept)", "Liking"), col.order = c("est", "se","ci",  "stat","df", "p"))

	BOLD Activity Value Network
Predictors	Estimates	CI	t	df	p
(Intercept)	-0.01	-0.09 – 0.08	-0.17	30.00	0.868
Liking	0.16	0.02 – 0.31	2.19	25.00	0.038
Random Effects
σ²	0.96
τ₀₀ _{subj_idx}	0.05
τ₁₁ _{subj_idx.cLiking}	0.01
ρ₀₁ _{subj_idx}	1.00
ICC	0.05
N _{subj_idx}	29
Observations	4126
Marginal R² / Conditional R²	0.001 / 0.047

Constituting regions of the value network

BWvMPFCLikMod <- lmer(BARTRA_vmPFC~ cLiking+(1|subj_idx), data = a1fMRI,
                                     REML=FALSE)

sv1_max <- svd(getME(BWvMPFCLikMod, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

BWStrLikMod <- lmer(BARTRA_NAcc~ cLiking+(1|subj_idx), data = a1fMRI,
                                   REML=FALSE)

sv1_max <- svd(getME(BWStrLikMod, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(BWvMPFCLikMod, BWStrLikMod,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("BOLD Activity vmPFC", "BOLD Activity vStr"), pred.labels= c("(Intercept)", "Liking"), col.order = c("est", "se","ci",  "stat","df", "p"))

	BOLD Activity vmPFC					BOLD Activity vStr
Predictors	Estimates	CI	t	df	p	Estimates	CI	t	df	p
(Intercept)	-0.57	-0.72 – -0.42	-7.43	30.00	<0.001	-0.42	-0.56 – -0.28	-5.89	30.00	<0.001
Liking	0.22	-0.03 – 0.48	1.70	3713.00	0.089	0.25	0.04 – 0.46	2.33	3901.00	0.020
Random Effects
σ²	3.29					2.21
τ₀₀	0.14 _{subj_idx}					0.13 _{subj_idx}
ICC	0.04					0.05
N	29 _{subj_idx}					29 _{subj_idx}
Observations	4126					4126
Marginal R² / Conditional R²	0.001 / 0.042					0.002 / 0.056

Is there an interaction with region?

# test for differences for Liking
BARTRASubLikMod0 <- lmer(BOLD~ cLiking*Region+(cLiking+Region|subj_idx), data = a2R,
                                     REML=FALSE)
sv1_max <- svd(getME(BARTRASubLikMod0, "Tlist")[[1]])
varience_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(BARTRASubLikMod0,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("BOLD Activity vmPFC vs vStr"), pred.labels= c("(Intercept)", "Liking", "Region (linear)", "Region: Liking"), col.order = c("est", "se","ci",  "stat","df", "p"))

	BOLD Activity vmPFC vs vStr
Predictors	Estimates	CI	t	df	p
(Intercept)	-0.49	-0.62 – -0.37	-7.83	30.00	<0.001
Liking	0.25	0.05 – 0.45	2.42	28.00	0.022
Region (linear)	0.11	-0.00 – 0.21	1.88	30.00	0.070
Region: Liking	0.01	-0.22 – 0.25	0.12	4639.00	0.901
Random Effects
σ²	2.75
τ₀₀ _{subj_idx}	0.10
τ₁₁ _{subj_idx.cLiking}	0.07
τ₁₁ _{subj_idx.Region.L}	0.07
ρ₀₁	0.42
	-0.08
ICC	0.05
N _{subj_idx}	29
Observations	8252
Marginal R² / Conditional R²	0.003 / 0.051

No, there is no interaction with region.

Follow up Striatal gradient for set appraisal

Figure S6

library(lmerTest)
BWFCN_StrCallInModLik <- lmer(FCN_Str~ cLiking+(cLiking|subj_idx), data = a1fMRI,
                                          REML=FALSE)

sv1_max <- svd(getME(BWFCN_StrCallInModLik, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

## make table part for liking figure
TEST <- as.data.frame(coef(summary(BWFCN_StrCallInModLik)))
TEST$effects <- rownames(TEST)
TEST$se <- TEST[,2]
TEST$effects<- revalue(TEST$effects, c("(Intercept)"="Intercept", "cLiking"= "Liking"))
TEST$effects<- ordered(TEST$effects, levels=c("Intercept", "Liking"))
TEST <- TEST[c(2), ]
TEST$effects <- factor(TEST$effects) 
TESTDCLik <- TEST

BWDMN_StrCallInModLik <- lmer(DMN_Str~ cLiking+(cLiking|subj_idx), data = a1fMRI,
                                          REML=FALSE)

sv1_max <- svd(getME(BWDMN_StrCallInModLik, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

TEST <- as.data.frame(coef(summary(BWDMN_StrCallInModLik)))
TEST$effects <- rownames(TEST)
TEST$se <- TEST[,2]
TEST$effects<- revalue(TEST$effects, c("(Intercept)"="Intercept", "cLiking"= "Liking"))
TEST$effects<- ordered(TEST$effects, levels=c("Intercept", "Liking"))
TEST <- TEST[c(2), ]
TEST$effects <- factor(TEST$effects) 
TESTVSSLik <- TEST

BWLN_StrCallInModLik <- lmer(LN_Str~ cLiking+(cLiking|subj_idx), data = a1fMRI,
                                         REML=FALSE)

sv1_max <- svd(getME(BWLN_StrCallInModLik, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

TEST <- as.data.frame(coef(summary(BWLN_StrCallInModLik)))
TEST$effects <- rownames(TEST)
TEST$se <- TEST[,2]
TEST$effects<- revalue(TEST$effects, c("(Intercept)"="Intercept", "cLiking"= "Liking"))
TEST$effects<- ordered(TEST$effects, levels=c("Intercept", "Liking"))
TEST <- TEST[c(2), ]
TEST$effects <- factor(TEST$effects) 
TESTVSILik <- TEST

detach("package:lmerTest", unload = TRUE)

## add region to table before merging
TESTDCLik$Region <- rep("Dorsal\n Caudate", length(TESTDCLik$Estimate))
TESTVSSLik$Region <- rep("Ventral\n Striatum\n (superior)", length(TESTVSSLik$Estimate))
TESTVSILik$Region <- rep("Ventral\n Striatum\n (inferior)", length(TESTVSILik$Estimate))
TEST2 <- rbind(TESTDCLik, TESTVSSLik, TESTVSILik)
TEST2$Region<- ordered(TEST2$Region, levels=c("Dorsal\n Caudate", "Ventral\n Striatum\n (superior)", "Ventral\n Striatum\n (inferior)"))

TEST2$isSig <- rep("*", length(TEST2$Estimate))
TEST2$isSig[TEST2[,5]>0.05] <- NA_real_ 
TEST2$isSig2 <- TEST2$isSig
TEST2$isSig2[TEST2[,5]<0.01]<- "**" 
TEST2$isSig2[TEST2[,5]<0.001]<- "***" 
EFFPLOTIA <- ggplot(data=TEST2, aes(x=Region, y=Estimate, fill=effects, group=effects, color=effects, shape=effects)) + geom_point(size=7)+theme_bw(14)+ scale_color_manual(values=c("#FF9900", "#70AD47", "#88d356"))+ylim(-0.01,0.6)+
  geom_errorbar(aes(max = Estimate + se, min = Estimate -se), width=0.1)+scale_shape_manual(values = c( 15, 16,18))+
  xlab("") + ylab("Regression Coefficient") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="none")+
  geom_hline(aes(yintercept=0), size=0.3,linetype=2, color="black")+ 
  geom_text(aes(label=isSig2, y=0.58), size=12,show.legend = F) 

EFFPLOTIA

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/BASBWfMRI/BARTRAROIEFFECTS_Liking_all_Regions_Str.pdf", width = 6, height = 4)
EFFPLOTIA
dev.off()

pdf 
  2

4) Supplemental Analyses

PCC

Table S4

BWBACallInModPCC <- lmer(POS_PCC~ isChooseBestTrial+cAvBid+ ChosenvAvRemV +cavGoalValue+ crelGoalvAvRem+cRT+(isChooseBestTrial+cRT|subj_idx), data = a1fMRI,
                                      REML=FALSE)

sv1_max <- svd(getME(BWBACallInModPCC, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(BWBACallInModPCC,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("BOLD Activity PCC"), pred.labels= c("(Intercept)", "Best - Worst Condition" ,"Overall Reward Value", "Relative Reward Value", "Overall Goal Value", "Relative Goal Value", "RT"), col.order = c("est", "se","ci",  "stat","df", "p"))

	BOLD Activity PCC
Predictors	Estimates	CI	t	df	p
(Intercept)	-0.76	-0.90 – -0.63	-11.02	31.00	<0.001
Best - Worst Condition	0.04	-0.13 – 0.22	0.51	41.00	0.616
Overall Reward Value	0.25	0.03 – 0.48	2.21	3766.00	0.027
Relative Reward Value	-0.09	-0.30 – 0.13	-0.77	4216.00	0.443
Overall Goal Value	0.12	-0.11 – 0.34	1.03	1829.00	0.301
Relative Goal Value	0.43	0.21 – 0.66	3.77	3999.00	<0.001
RT	-0.16	-0.32 – -0.01	-2.08	31.00	0.046
Random Effects
σ²	2.81
τ₀₀ _{subj_idx}	0.12
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.12
τ₁₁ _{subj_idx.cRT}	0.04
ρ₀₁	0.16
	-0.01
ICC	0.05
N _{subj_idx}	30
Observations	4270
Marginal R² / Conditional R²	0.008 / 0.060

rACC vs mOFC

Is there an interaction with region within vmPFC?

Table S6

a2RVMPFC <- rbind(a1fMRI, a1fMRI)

a2RVMPFC$Region <- rep(c("rACC", "OFC"), each=length(a1fMRI$subj_idx))
a2RVMPFC$Region <- as.factor(a2RVMPFC$Region)
a2RVMPFC$Region <- ordered(a2RVMPFC$Region, levels= c("rACC", "OFC"))
contrasts(a2RVMPFC$Region) <- contr.sdif(2)
a2RVMPFC$BOLD <- a2RVMPFC$sPOS_rACC
a2RVMPFC$BOLD[a2RVMPFC$Region=="OFC"] <- a2RVMPFC$sHUvC_mOFC[a2RVMPFC$Region=="OFC"] 
contrasts(a2RVMPFC$isChooseBestTrial) <- contr.sdif(2)
library(lmerTest)
# no significant interactions with Region
VMPFCSubMod0 <- lmer(BOLD~ (isChooseBestTrial+cAvBid+ ChosenvAvRemV +cavGoalValue+ crelGoalvAvRem)*Region+cRT+(isChooseBestTrial+cRT+Region|subj_idx), data = a2RVMPFC,
                                     REML=FALSE)

sv1_max <- svd(getME(VMPFCSubMod0, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

xxx <- anova(VMPFCSubMod0)
row.names(xxx) <- c("Best - Worst Condition","Overall Reward Value", "Relative Reward Value", "Overall Goal Value", "Relative Goal Value", "Region", "RT","Region: Best - Worst Condition", "Region: Overall Reward Value", "Region: Relative Reward Value", "Region: Overall Goal Value", "Region: Relative Goal Value" )
attributes(xxx)$heading <- NULL
xxx[, c(1:5)] <- round(xxx[, c(1:5)],2)
xxx[, 6] <- round(xxx[, 6],3)
htmlTable(xxx)

	Sum Sq	Mean Sq	NumDF	DenDF	F value	Pr(>F)
Best - Worst Condition	1.6	1.6	1	226.88	1.67	0.197
Overall Reward Value	5.38	5.38	1	6676.4	5.63	0.018
Relative Reward Value	0.06	0.06	1	8475.2	0.07	0.796
Overall Goal Value	2.58	2.58	1	5966.15	2.7	0.101
Relative Goal Value	6.61	6.61	1	8405.74	6.91	0.009
Region	0.02	0.02	1	29.39	0.02	0.89
RT	0.23	0.23	1	28.54	0.24	0.628
Region: Best - Worst Condition	0.87	0.87	1	8440.16	0.91	0.341
Region: Overall Reward Value	0.21	0.21	1	2585.16	0.22	0.641
Region: Relative Reward Value	0.91	0.91	1	8455.28	0.95	0.33
Region: Overall Goal Value	0.05	0.05	1	8445.82	0.06	0.811
Region: Relative Goal Value	1.29	1.29	1	8358.15	1.34	0.246

detach("package:lmerTest", unload = TRUE)

No, there is no significant interaction.

Table S7

BWBACallInModrACC <- lmer(POS_rACC~ isChooseBestTrial+cAvBid+ ChosenvAvRemV +cavGoalValue+ crelGoalvAvRem+cRT+(cRT|subj_idx), data = a1fMRI,
                                        REML=FALSE)

sv1_max <- svd(getME(BWBACallInModrACC, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

BWBACallInModmOFC <- lmer(HUvC_mOFC~ isChooseBestTrial+cAvBid+ ChosenvAvRemV +cavGoalValue+ crelGoalvAvRem+cRT +(cRT|subj_idx), data = a1fMRI,
                                        REML=FALSE)

sv1_max <- svd(getME(BWBACallInModmOFC, "Tlist")[[1]])
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(BWBACallInModrACC,BWBACallInModmOFC,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("BOLD Activity rACC", "BOLD Activity mOFC"), pred.labels= c("(Intercept)", "Best - Worst Condition" ,"Overall Reward Value", "Relative Reward Value", "Overall Goal Value", "Relative Goal Value", "RT"), col.order = c("est", "se","ci",  "stat","df", "p"))

	BOLD Activity rACC					BOLD Activity mOFC
Predictors	Estimates	CI	t	df	p	Estimates	CI	t	df	p
(Intercept)	-0.77	-0.87 – -0.67	-14.85	31.00	<0.001	-0.23	-0.34 – -0.12	-3.95	31.00	<0.001
Best - Worst Condition	-0.02	-0.13 – 0.09	-0.32	4087.00	0.751	-0.08	-0.20 – 0.03	-1.44	4154.00	0.150
Overall Reward Value	0.22	0.02 – 0.43	2.15	3253.00	0.032	0.12	-0.09 – 0.33	1.16	3562.00	0.247
Relative Reward Value	-0.06	-0.26 – 0.14	-0.61	4252.00	0.544	0.10	-0.11 – 0.31	0.95	4249.00	0.341
Overall Goal Value	0.02	-0.17 – 0.21	0.22	4216.00	0.830	0.22	0.03 – 0.42	2.24	4240.00	0.025
Relative Goal Value	0.18	-0.03 – 0.39	1.70	4203.00	0.089	0.24	0.03 – 0.45	2.19	4249.00	0.028
RT	-0.36	-0.53 – -0.18	-4.04	32.00	<0.001	0.30	0.10 – 0.49	2.98	33.00	0.005
Random Effects
σ²	2.37					2.50
τ₀₀	0.06 _{subj_idx}					0.08 _{subj_idx}
τ₁₁	0.11 _{subj_idx.cRT}					0.17 _{subj_idx.cRT}
ρ₀₁	0.17 _{subj_idx}					0.04 _{subj_idx}
ICC	0.03					0.04
N	30 _{subj_idx}					30 _{subj_idx}
Observations	4270					4270
Marginal R² / Conditional R²	0.014 / 0.047					0.009 / 0.052

Study S0: choose best only

the original effects: RT in choose best

Data Prep

a0$subj_idx <- as.factor(a0$subj_idx)
a0$cgoalvAvRemBid <- scale(a0$MaxvAvRemBid/10, center=TRUE, scale = FALSE)
a0$cAvBid <- scale(a0$AvBid/10, center=TRUE, scale = FALSE)
a0$logRT <- log(a0$rt)
a0$ChosenvAvRemV <- scale(a0$ChosenvAvRemV/10, scale=FALSE, center=TRUE)

accuracy

VD > 0 only

Baccmod0 <- glmer(response~ cgoalvAvRemBid+cAvBid+(1|subj_idx), data = a0[!a0$goalvNext==0,], 
                                 family = binomial)

sv1_max <- svd(getME(Baccmod0, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(Baccmod0, transform= NULL ,show.stat = TRUE, string.stat="z",   dv.labels = c("P(Max Chosen)"), pred.labels= c("(Intercept)", "Value Difference", "Overall Value"), col.order = c("est", "se","ci",  "stat","df", "p"))

	P(Max Chosen)
Predictors	Log-Odds	CI	z	p
(Intercept)	0.09	-0.17 – 0.36	0.68	0.496
Value Difference	2.98	2.29 – 3.67	8.48	<0.001
Overall Value	-0.63	-1.38 – 0.11	-1.66	0.097
Random Effects
σ²	3.29
τ₀₀ _{subj_idx}	0.22
ICC	0.06
N _{subj_idx}	19
Observations	1130
Marginal R² / Conditional R²	0.126 / 0.180

eff_df <- Effect(c("cgoalvAvRemBid"), Baccmod0,xlevels=list(cgoalvAvRemBid=seq(min(a0$cgoalvAvRemBid), max(a0$cgoalvAvRemBid), 0.1)))
contmain <- as.data.frame(eff_df)
contmain$cgoalvAvRemBid <- (contmain$cgoalvAvRemBid - min(contmain$cgoalvAvRemBid))*10
pbacc <- ggplot(data=contmain, aes(x=cgoalvAvRemBid, y=fit)) + geom_line(color="#000000")+theme_bw(12)+ geom_ribbon(data=contmain, aes(x=cgoalvAvRemBid, max = upper, min = lower),alpha=0.3,  inherit.aes = FALSE, show.legend=FALSE)+ #ggtitle("Choose worst")+theme(plot.title = element_text(hjust = 0.5))+ 
  scale_y_continuous(limits=c(0,1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+ggtitle("choose best")+
  xlab("Value difference") + ylab("Pr(max chosen)") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")#+ theme(axis.title.x=element_text(hjust= -5))
pbacc

RT

BRTmod0 <- lmer(logRT~ cgoalvAvRemBid+cAvBid+(cAvBid|subj_idx), data = a0,
                               REML=FALSE)

sv1_max <- svd(getME(BRTmod0, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(BRTmod0,p.val= "kr", show.df = T, transform= NULL ,show.stat = TRUE, string.stat="t",   dv.labels = c("log RT"),pred.labels= c("(Intercept)", "Value Difference", "Overall Value"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.39	1.28 – 1.50	24.71	20.00	<0.001
Value Difference	-0.67	-0.75 – -0.59	-15.70	2264.00	<0.001
Overall Value	-0.33	-0.45 – -0.20	-5.02	19.00	<0.001
Random Effects
σ²	0.23
τ₀₀ _{subj_idx}	0.05
τ₁₁ _{subj_idx.cAvBid}	0.04
ρ₀₁ _{subj_idx}	0.54
ICC	0.20
N _{subj_idx}	19
Observations	2280
Marginal R² / Conditional R²	0.098 / 0.282

eff_df <- Effect(c("cAvBid"), BRTmod0)
contmain <- as.data.frame(eff_df)
contmain$cAvBid <- (contmain$cAvBid +0.5)*10
pASVB <- ggplot(data=contmain, aes(x=cAvBid, y=fit)) + geom_line(color="#000000")+theme_bw(12)+ geom_ribbon(data=contmain, aes(x=cAvBid, max = upper, min = lower),alpha=0.3, inherit.aes = FALSE)+ ylim(0.7, 2)+theme(plot.title = element_text(hjust = 0.5))+ #ggtitle("Choose best")+ 
  scale_y_continuous(limits=c(0.7, 2), expand=c(0, 0))+ scale_x_continuous(expand=c(0, 0))+xlab("Overall Value ") + ylab("log RT") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")


eff_df <- Effect(c("cgoalvAvRemBid"), BRTmod0, xlevels=list(cgoalvAvRemBid=seq(min(a0$cgoalvAvRemBid), max(a0$cgoalvAvRemBid), 0.1) ))
contmain <- as.data.frame(eff_df)
#contmain$cgoalvAvRemBid <- contmain$cgoalvAvRemBid - min(contmain$cgoalvAvRemBid)

pVDB <- ggplot(data=contmain, aes(x=cgoalvAvRemBid, y=fit)) + geom_line(color="#000000")+theme_bw(12)+ geom_ribbon(data=contmain, aes(x=cgoalvAvRemBid, max = upper, min = lower),alpha=0.3, inherit.aes = FALSE)+ ylim(0.7, 2)+ ggtitle("")+
  scale_y_continuous(limits=c(0.7, 2), expand=c(0, 0))+ scale_x_continuous( expand=c(0, 0))+xlab("Relative Value") + ylab("log RT") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")
multiplot(pVDB, pASVB, cols=1)

Study S1: choose worst only

Data Prep

# 1: center covariates
a2$cgoalvAvRemBid <- scale(a2$goalvAvRemBid/10, center=TRUE, scale = FALSE)
a2$cAvBid <- scale(a2$AvBid/10, center=TRUE, scale = FALSE)
a2$ChosenvAvRemV <- scale(a2$ChosenvAvRemV/10, scale=FALSE, center=TRUE)   
# 2: make categorical predictors factors and set contrasts
a2$subj_idx <- as.factor(a2$subj_idx)

accuracy

Waccmod0 <- glmer(response~ cgoalvAvRemBid+cAvBid+(1|subj_idx), data = a2[!a2$goalvNext==0,], 
                                 family = binomial)

sv1_max <- svd(getME(Waccmod0, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(Waccmod0, transform= NULL ,show.stat = TRUE, string.stat="z",   dv.labels = c("P(Min Chosen)"), pred.labels= c("(Intercept)", "Value Difference", "Overall Value"), col.order = c("est", "se","ci",  "stat","df", "p"))

	P(Min Chosen)
Predictors	Log-Odds	CI	z	p
(Intercept)	-0.22	-0.40 – -0.03	-2.31	0.021
Value Difference	3.66	2.92 – 4.39	9.72	<0.001
Overall Value	0.05	-0.72 – 0.81	0.12	0.903
Random Effects
σ²	3.29
τ₀₀ _{subj_idx}	0.04
ICC	0.01
N _{subj_idx}	18
Observations	1097
Marginal R² / Conditional R²	0.130 / 0.140

eff_df <- Effect(c("cgoalvAvRemBid"), Waccmod0,xlevels=list(cgoalvAvRemBid=seq(min(a2$cgoalvAvRemBid[!a2$goalvNext==0]), max(a2$cgoalvAvRemBid[!a2$goalvNext==0]), 0.1)))
contmain <- as.data.frame(eff_df)
contmain$cgoalvAvRemBid <- (contmain$cgoalvAvRemBid - min(contmain$cgoalvAvRemBid))*10
pwacc <- ggplot(data=contmain, aes(x=cgoalvAvRemBid, y=fit)) + geom_line(color="#CC0000")+theme_bw(12)+ geom_ribbon(data=contmain, aes(x=cgoalvAvRemBid, max = upper, min = lower), fill="#CC0000",alpha=0.3,  inherit.aes = FALSE, show.legend=FALSE)+ #ggtitle("Choose worst")+theme(plot.title = element_text(hjust = 0.5))+ 
  scale_y_continuous(limits=c(0,1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+ ggtitle("choose worst")+
  xlab("Value difference") + ylab("Pr(min chosen)") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")#+ theme(axis.title.x=element_text(hjust= -5))

multiplot(pbacc, pwacc, cols = 2)

pdf("/Users/Romy/Dropbox (Brown)/ShenhavLab/Experiments/bas/Analysis/HDDM/best_worst_between_accuracy.pdf", width = 7, height = 4)#, units = 'cm', res = 200, compression = 'lzw'
multiplot(pbacc, pwacc, cols = 2)
dev.off()

pdf 
  2

rt

WRTmod0 <- lmer(logRT~ cgoalvAvRemBid+cAvBid+(cAvBid|subj_idx), data = a2,
                               REML=FALSE)

sv1_max <- svd(getME(WRTmod0, "Tlist")[[1]]) 
variance_explainned <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)

tab_model(WRTmod0,p.val= "kr", show.df = T, transform= NULL ,show.stat = TRUE, string.stat="t",   dv.labels = c("log RT"),pred.labels= c("(Intercept)", "Value Difference", "Overall Value"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.50	1.34 – 1.67	18.20	19.00	<0.001
Value Difference	-0.57	-0.67 – -0.47	-11.09	2122.00	<0.001
Overall Value	0.52	0.32 – 0.72	4.99	19.00	<0.001
Random Effects
σ²	0.26
τ₀₀ _{subj_idx}	0.11
τ₁₁ _{subj_idx.cAvBid}	0.13
ρ₀₁ _{subj_idx}	-0.20
ICC	0.32
N _{subj_idx}	18
Observations	2148
Marginal R² / Conditional R²	0.066 / 0.361

eff_df <- Effect(c("cAvBid"), WRTmod0)
contmain <- as.data.frame(eff_df)
contmain$cAvBid <- (contmain$cAvBid +0.5)*10
pASVW <- ggplot(data=contmain, aes(x=cAvBid, y=fit)) + geom_line(color="#CC0000")+theme_bw(12)+ geom_ribbon(data=contmain, aes(x=cAvBid, max = upper, min = lower), fill="#CC0000",alpha=0.3,  inherit.aes = FALSE, show.legend=FALSE)+ ylim(0.7, 2)+theme(plot.title = element_text(hjust = 0.5))+ #ggtitle("Choose worst")+
  scale_y_continuous(limits=c(0.7, 2), expand=c(0, 0))+ scale_x_continuous(expand=c(0, 0))+xlab(" ") + ylab(" ") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")#+ theme(axis.title.x=element_text(hjust= -5))


eff_df <- Effect(c("cgoalvAvRemBid"), WRTmod0, xlevels=list(cgoalvAvRemBid=seq(min(a2$cgoalvAvRemBid), max(a2$cgoalvAvRemBid), 0.1) ))
contmain <- as.data.frame(eff_df)
#contmain$cgoalvAvRemBid <- contmain$cgoalvAvRemBid - min(contmain$cgoalvAvRemBid)
pVDW <- ggplot(data=contmain, aes(x=cgoalvAvRemBid, y=fit)) + geom_line(color="#CC0000")+theme_bw(12)+ geom_ribbon(data=contmain, aes(x=cgoalvAvRemBid, max = upper, min = lower), fill="#CC0000",alpha=0.3, inherit.aes = FALSE, show.legend=FALSE)+ylim(0.7, 2)+ ggtitle("")+
  scale_y_continuous(limits=c(0.7, 2), expand=c(0, 0))+ scale_x_continuous(expand=c(0, 0))+xlab(" ") + ylab(" ") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="bottom")#+ theme(axis.title.x=element_text(hjust= -5))

multiplot(pVDB, pASVB, pVDW, pASVW,  cols=2)

pdf("/Users/Romy/Dropbox (Brown)/ShenhavLab/Experiments/bas/Analysis/HDDM/best_worst_between.pdf", width = 7, height = 8)#, units = 'cm', res = 200, compression = 'lzw'
multiplot(pVDB, pASVB, pVDW, pASVW,  cols=2)
dev.off()

pdf 
  2

Study S2: best vs worst within, incentivized

data prep

# 1: center covariates
a118$cgoalvAvRemBid <- scale(a118$goalvAvRemBid/10, center=TRUE, scale = FALSE)                  ## unsigned value difference (max/min vs average remaining)
a118$cgoalvNext <- scale(a118$goalvNext, center=TRUE, scale = FALSE)
a118$cAvBid <- scale(a118$AvBid/10, center=TRUE, scale = FALSE)
a118$cRT <- scale(a118$logRT, center=TRUE, scale=FALSE)
#a118$cChosenvAvRemV <- scale(a118$ChosenvAvRemV, center= TRUE, scale=FALSE)
a118$Liking <- scale(a118$tmpAllLrateAllTrials, scale=FALSE, center=TRUE)
a118$avGoalValue <- a118$AvBid                                                                ## the value of the set relative to the goal
a118$avGoalValue[a118$isChooseBestTrial==0] <- (a118$AvBid [a118$isChooseBestTrial==0]*-1) +10
a118$cavGoalValue <- scale(a118$avGoalValue/10, scale=FALSE, center= TRUE)
a118$goalChosenV <- a118$ChosenV                                                              ## the value of the chosen Item relativ to goal
a118$goalChosenV[a118$isChooseBestTrial==0] <- (a118$ChosenV[a118$isChooseBestTrial==0]*-1) +10
a118$cgoalChosenV <- scale(a118$goalChosenV, scale=FALSE, center= TRUE)
a118$cChosenV <- scale(a118$ChosenV, scale=FALSE, center=TRUE)
a118$avUnchV <- a118$SumUnchV/3                                                               ## average of the unchosen values
a118$cavUnchV <- scale(a118$avUnchV, scale=FALSE, center=TRUE)
a118$avgoalUnchV <- a118$avUnchV                                                              ## average of unchosen values relative to goal
a118$avgoalUnchV[a118$isChooseBestTrial==0] <- (a118$avUnchV[a118$isChooseBestTrial==0]*-1) +10
a118$cavgoalUnchV <- scale(a118$avgoalUnchV, scale=FALSE, center=TRUE)
a118$relGoalvAvRem <- a118$goalChosenV - a118$avgoalUnchV                                       ## goal related chosen vs unchosen (negative Values are suboptimal choices) 
a118$crelGoalvAvRem <- scale(a118$relGoalvAvRem/10, scale=FALSE, center=TRUE)
a118$goalVal <- a118$MaxBid
a118$goalVal[a118$isChooseBestTrial==0] <- a118$MinBid[a118$isChooseBestTrial==0]
a118$cgoalVal <- scale(a118$goalVal, scale=FALSE, center=TRUE)                               ## goal Value (best matching value in set relative to goal)
a118$rewSpgoalvAvRem <- a118$goalVal - a118$avUnchV                                            ## goal v av rem in bid space
a118$crewSpgoalvAvRem <- scale(a118$rewSpgoalvAvRem/10, scale=FALSE, center=TRUE) 
a118$Prgoalcons <- a118$goalVal
a118$Prgoalcons[a118$isChooseBestTrial==0] <- 10 - a118$MinBid[a118$isChooseBestTrial==0]        ## this is a terrible variable name, but that's the revalued goal value
a118$cChosenvAvRemV <- scale(a118$ChosenvAvRemV/10, scale=FALSE, center=TRUE)   
a118$chosenvsnext <- a118$ChosenV - a118$allTrProdBid2

# 2: make categorical predictors factors and set contrasts
a118$subj_idx <- as.factor(a118$subj_idx)
a118$isChooseBestTrial <- as.factor(a118$isChooseBestTrial)
contrasts(a118$isChooseBestTrial) <- contr.sdif(2)
a118$block <- rep(1, length(a118$subj_idx))
a118$block[a118$Trial > 60] <- 2
a118$block <- as.factor(a118$block)
contrasts(a118$block) <- contr.sdif(2)
a118$CondType <- as.factor(a118$CondType)
a118$CondType <- revalue(a118$CondType, c("1"="mixed", "2"="low", "3"="medium", "4"="high"))


## exclude trials with RT longer than 30s
as <- a118
a118 <- a118[a118$rt<30 & !a118$subj_idx==1201,]

accuracy

BWaccmod0t1 <- glmer(response~ cgoalvAvRemBid*isChooseBestTrial+cAvBid+(1|subj_idx), data = a118[!a118$goalvNext==0,], #
                                    family = binomial)

sv1_max <- svd(getME(BWaccmod0, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  


tab_model(BWaccmod0t1, transform= NULL ,show.stat = TRUE, string.stat="z",   dv.labels = c("P(Goal Chosen)"), pred.labels= c("(Intercept)", "Value Difference","Best - Worst Conditon", "Overall Value", "Value Difference:Best - Worst Conditon"), col.order = c("est", "se","ci",  "stat","df", "p"))

	P(Goal Chosen)
Predictors	Log-Odds	CI	z	p
(Intercept)	-0.16	-0.40 – 0.08	-1.33	0.183
Value Difference	2.53	1.77 – 3.30	6.50	<0.001
Best - Worst Conditon	0.27	-0.11 – 0.65	1.40	0.161
Overall Value	0.50	-0.29 – 1.29	1.24	0.214
Value Difference:Best - Worst Conditon	-0.84	-2.35 – 0.67	-1.09	0.274
Random Effects
σ²	3.29
τ₀₀ _{subj_idx}	0.07
ICC	0.02
N _{subj_idx}	14
Observations	799
Marginal R² / Conditional R²	0.088 / 0.107

no significant interactions, so excluding interaction term

BWaccmod0 <- glmer(response~ cgoalvAvRemBid+isChooseBestTrial+cAvBid+(1|subj_idx), data = a118[!a118$goalvNext==0,], #
                                  family = binomial)

sv1_max <- svd(getME(BWaccmod0, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(BWaccmod0, transform= NULL ,show.stat = TRUE, string.stat="z",   dv.labels = c("P(Goal Chosen)"), pred.labels= c("(Intercept)", "Value Difference","Best - Worst Conditon", "Overall Value"), col.order = c("est", "se","ci",  "stat","df", "p"))

	P(Goal Chosen)
Predictors	Log-Odds	CI	z	p
(Intercept)	-0.16	-0.40 – 0.08	-1.34	0.182
Value Difference	2.54	1.78 – 3.30	6.55	<0.001
Best - Worst Conditon	0.14	-0.16 – 0.44	0.92	0.357
Overall Value	0.57	-0.21 – 1.35	1.44	0.150
Random Effects
σ²	3.29
τ₀₀ _{subj_idx}	0.07
ICC	0.02
N _{subj_idx}	14
Observations	799
Marginal R² / Conditional R²	0.085 / 0.104

eff_df <- Effect(c("cgoalvAvRemBid", "isChooseBestTrial"),BWaccmod0, xlevels=list(cgoalvAvRemBid=seq(min(a118$cgoalvAvRemBid), max(a118$cgoalvAvRemBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))
IA$cgoalvAvRemBid <- (IA$cgoalvAvRemBid - min(IA$cgoalvAvRemBid))*10

## with CI instead of se
pBWacc <- ggplot(data=IA, aes(x=cgoalvAvRemBid, y=fit , color=isChooseBestTrial )) + geom_line()+scale_colour_manual(name="Goal",values=c("#CC0000","#000000")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cgoalvAvRemBid, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous( expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("value difference") + ylab("Pr(goal chosen)") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.justification=c(1,0), legend.position=c(0.95,0.05))

pBWacc

pdf("/Users/Romy/Dropbox (Brown)/ShenhavLab/Experiments/bas/Analysis/HDDM/best_worst_within_main_accuracy_w_choice_condition_as_color118.pdf", width = 4, height = 4)#, units = 'cm', res = 200, compression = 'lzw'
pBWacc 
dev.off()

pdf 
  2

S1BWaccmod0G <- glmer(response~ cgoalvAvRemBid*isChooseBestTrial+cavGoalValue+(isChooseBestTrial|subj_idx), data = a118[!a118$goalvNext==0,], 
                                  family = binomial) 

sv1_max <- svd(getME(S1BWaccmod0G, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(S1BWaccmod0G, transform= NULL ,show.stat = TRUE, string.stat="z",   dv.labels = c("P(Goal Chosen)"), pred.labels= c("(Intercept)", "Value Difference","Best - Worst Conditon", "Overall Goal Value", "Value Difference:Best - Worst Conditon"), col.order = c("est", "se","ci",  "stat","df", "p"))

	P(Goal Chosen)
Predictors	Log-Odds	CI	z	p
(Intercept)	-0.14	-0.32 – 0.05	-1.41	0.158
Value Difference	2.40	1.66 – 3.15	6.32	<0.001
Best - Worst Conditon	0.24	-0.16 – 0.64	1.16	0.247
Overall Goal Value	-0.44	-1.17 – 0.28	-1.20	0.230
Value Difference:Best - Worst Conditon	-0.93	-2.41 – 0.55	-1.23	0.219
Random Effects
σ²	3.29
τ₀₀ _{subj_idx}	0.00
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.03
ρ₀₁ _{subj_idx}
N _{subj_idx}	14
Observations	799
Marginal R² / Conditional R²	0.084 / NA

rt

Basic checks

meanRTbySubS1 <- ddply(a118, .(subj_idx, isChooseBestTrial), summarise,
                       mrt    = median(logRT, na.rm=TRUE))


ggplot(meanRTbySubS1, aes(isChooseBestTrial,mrt,  color= isChooseBestTrial)) + geom_violin()+geom_boxplot() + theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))

ggplot(a118, aes(rt, group =isChooseBestTrial, color=isChooseBestTrial, fill=isChooseBestTrial)) + geom_density(alpha=0.5) + theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))

analyses separately for best and worst before putting all in one model

BWRTmod1bBEST <- lmer(logRT~ cgoalvAvRemBid+cAvBid+(1|subj_idx), data = a118[a118$isChooseBestTrial==1,],
                                     REML=FALSE)

sv1_max <- svd(getME(BWRTmod1bBEST, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  


BWRTmod1bWORST <- lmer(logRT~ cgoalvAvRemBid+cAvBid+(1|subj_idx), data = a118[a118$isChooseBestTrial==0,],
                                      REML=FALSE)

sv1_max <- svd(getME(BWRTmod1bWORST, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(BWRTmod1bBEST,BWRTmod1bWORST, p.val= "kr", show.df = T, transform= NULL ,show.stat = TRUE, string.stat="t",   dv.labels = c("log RT Best", "log RT worst"),pred.labels= c("(Intercept)", "Value Difference", "Overall Value"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT Best					log RT worst
Predictors	Estimates	CI	t	df	p	Estimates	CI	t	df	p
(Intercept)	1.74	1.57 – 1.91	19.90	15.00	<0.001	1.66	1.45 – 1.87	15.44	15.00	<0.001
Value Difference	-0.65	-0.80 – -0.50	-8.51	814.00	<0.001	-0.57	-0.70 – -0.44	-8.45	811.00	<0.001
Overall Value	-0.22	-0.36 – -0.08	-3.14	825.00	0.002	0.20	0.06 – 0.33	2.84	820.00	0.005
Random Effects
σ²	0.23					0.20
τ₀₀	0.09 _{subj_idx}					0.15 _{subj_idx}
ICC	0.29					0.43
N	14 _{subj_idx}					14 _{subj_idx}
Observations	823					821
Marginal R² / Conditional R²	0.072 / 0.340					0.053 / 0.459

Standard analysis without interaction of choice goal and overall value

# all data no interaction
BWRTmod0 <- lmer(logRT~ cgoalvAvRemBid+cAvBid+isChooseBestTrial+(isChooseBestTrial|subj_idx), data = a118,
                                REML=FALSE)

sv1_max <- svd(getME(BWRTmod0, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(BWRTmod0,p.val= "kr", show.df = T, transform= NULL ,show.stat = TRUE, string.stat="t",   dv.labels = c("log RT"),pred.labels= c("(Intercept)", "Value Difference", "Overall Value", "Best - Worst Condition"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.70	1.52 – 1.88	18.70	15.00	<0.001
Value Difference	-0.60	-0.70 – -0.50	-11.69	1629.00	<0.001
Overall Value	-0.02	-0.12 – 0.08	-0.42	1634.00	0.677
Best - Worst Condition	0.08	-0.07 – 0.22	1.07	15.00	0.300
Random Effects
σ²	0.22
τ₀₀ _{subj_idx}	0.11
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.06
ρ₀₁ _{subj_idx}	-0.34
ICC	0.36
N _{subj_idx}	14
Observations	1644
Marginal R² / Conditional R²	0.059 / 0.396

all data with interaction

Table S2

BWRTmod1b <- lmer(logRT~ cgoalvAvRemBid+cAvBid*isChooseBestTrial+(isChooseBestTrial|subj_idx), data = a118,
                                 REML=FALSE)

sv1_max <- svd(getME(BWRTmod1b, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(BWRTmod1b,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("log RT"), pred.labels= c("(Intercept)", "Value Difference", "Overall Value", "Best - Worst Condition", "Overall Value: Best - Worst"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.70	1.52 – 1.88	18.77	15.00	<0.001
Value Difference	-0.61	-0.71 – -0.51	-12.05	1630.00	<0.001
Overall Value	-0.01	-0.11 – 0.09	-0.20	1635.00	0.838
Best - Worst Condition	0.08	-0.07 – 0.22	1.00	15.00	0.335
Overall Value: Best - Worst	-0.41	-0.61 – -0.22	-4.20	1489.00	<0.001
Random Effects
σ²	0.21
τ₀₀ _{subj_idx}	0.10
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.07
ρ₀₁ _{subj_idx}	-0.30
ICC	0.36
N _{subj_idx}	14
Observations	1644
Marginal R² / Conditional R²	0.066 / 0.404

nested model to get by choice goal estimates of overall value

BWRTmod1bn <- lmer(logRT~ isChooseBestTrial/(cAvBid)+cgoalvAvRemBid+(isChooseBestTrial|subj_idx), data = a118,
                                  REML=FALSE)

sv1_max <- svd(getME(BWRTmod1bn, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  


tab_model(BWRTmod1bn,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("log RT"), pred.labels= c("(Intercept)", "Best - Worst Condition", "Value Difference", "Worst: Overall Value", "Best: Overall Value"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.70	1.52 – 1.88	18.77	15.00	<0.001
Best - Worst Condition	0.08	-0.07 – 0.22	1.00	15.00	0.335
Value Difference	-0.61	-0.71 – -0.51	-12.05	1630.00	<0.001
Worst: Overall Value	0.20	0.06 – 0.34	2.74	1614.00	0.006
Best: Overall Value	-0.22	-0.35 – -0.08	-3.18	1583.00	0.002
Random Effects
σ²	0.21
τ₀₀ _{subj_idx}	0.10
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.07
ρ₀₁ _{subj_idx}	-0.30
ICC	0.36
N _{subj_idx}	14
Observations	1644
Marginal R² / Conditional R²	0.066 / 0.404

eff_df <- Effect(c("cAvBid", "isChooseBestTrial"), BWRTmod1b)
IA <- as.data.frame(eff_df)
IA$cAvBid <- (IA$cAvBid - min(IA$cAvBid))*10
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
pASVBW <- ggplot(data=IA, aes(x=cAvBid, y=fit , color=isChooseBestTrial )) + geom_line()+scale_colour_manual(name="Goal",values=c("#CC0000","#000000")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cAvBid, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous(limits=c(0.8, 2.1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Overall value") + ylab(" ") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+ theme(legend.justification=c(1,1), legend.position=c(0.95,0.25) )


eff_df <- Effect(c("cgoalvAvRemBid", "isChooseBestTrial"), BWRTmod1b, xlevels=list(cgoalvAvRemBid=seq(min(a118$cgoalvAvRemBid), max(a118$cgoalvAvRemBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))
IA$cgoalvAvRemBid <- (IA$cgoalvAvRemBid - min(IA$cgoalvAvRemBid))*10
## with CI instead of se
pVDBW <- ggplot(data=IA, aes(x=cgoalvAvRemBid, y=fit , color=isChooseBestTrial )) + geom_line()+scale_colour_manual(name="Goal",values=c("#CC0000","#000000")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cgoalvAvRemBid, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous(limits=c(0.8, 2.1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Relative Value") + ylab("log RT") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+ theme(legend.position="none" )

multiplot(pVDBW, pASVBW,  cols=2)

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/HDDM/best_worst_within_allRS118.pdf", width = 8, height = 4)#, units = 'cm', res = 200, compression = 'lzw'
multiplot(pVDBW, pASVBW,  cols=2)
dev.off()

pdf 
  2

goal value instead of reward value

BWRTmod1g <- lmer(logRT~ cgoalvAvRemBid+cavGoalValue+isChooseBestTrial+(isChooseBestTrial|subj_idx), data = a118,
                                 REML=FALSE)

sv1_max <- svd(getME(BWRTmod1g, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  


tab_model(BWRTmod1g,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("log RT"), pred.labels= c("(Intercept)", "Value Difference", "Overall Goal Value", "Best - Worst Condition"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.70	1.52 – 1.88	18.79	15.00	<0.001
Value Difference	-0.61	-0.71 – -0.52	-12.13	1628.00	<0.001
Overall Goal Value	-0.21	-0.30 – -0.11	-4.21	1485.00	<0.001
Best - Worst Condition	0.05	-0.10 – 0.20	0.69	15.00	0.499
Random Effects
σ²	0.21
τ₀₀ _{subj_idx}	0.10
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.07
ρ₀₁ _{subj_idx}	-0.30
ICC	0.36
N _{subj_idx}	14
Observations	1644
Marginal R² / Conditional R²	0.066 / 0.404

eff_df <- Effect(c("cavGoalValue", "isChooseBestTrial"), BWRTmod1g)
IA <- as.data.frame(eff_df)
IA$cavGoalValue <-(IA$cavGoalValue - min(IA$cavGoalValue))*10
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
pASVBW <- ggplot(data=IA, aes(x=cavGoalValue, y=fit , color=isChooseBestTrial )) + geom_line()+scale_colour_manual(name="Goal",values=c("#CC0000","#000000")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cavGoalValue, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous(limits=c(0.8, 2.1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Overall value") + ylab(" ") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.justification=c(1,1), legend.position=c(0.95,0.25))

eff_df <- Effect(c("cgoalvAvRemBid", "isChooseBestTrial"), BWRTmod1g, xlevels=list(cgoalvAvRemBid=seq(min(a118$cgoalvAvRemBid), max(a118$cgoalvAvRemBid), 0.1)))
IA <- as.data.frame(eff_df)
IA$cgoalvAvRemBid <- (IA$cgoalvAvRemBid - min(IA$cgoalvAvRemBid))*10
IA$isChooseBestTrial<- revalue(IA$isChooseBestTrial, c("0"="worst", "1"="best"))

## with CI instead of se
pVDBW <- ggplot(data=IA, aes(x=cgoalvAvRemBid, y=fit , color=isChooseBestTrial )) + geom_line()+scale_colour_manual(name="Goal",values=c("#CC0000","#000000")) +theme_bw(12)+ geom_ribbon(data=IA, aes(x=cgoalvAvRemBid, max = upper, min = lower, fill = isChooseBestTrial),alpha=0.3, inherit.aes = FALSE)+
  scale_y_continuous(limits=c(0.8, 2.1), expand=c(0, 0))+ scale_x_continuous(limits=c(0, 10), expand=c(0, 0))+scale_fill_manual(name="Goal",values=c("#CC0000","#000000")) + xlab("Relative Value") + ylab("log RT") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position="none" )
multiplot( pVDBW,pASVBW, cols=2)

pdf("~/Dropbox (Brown)/ShenhavLab/experiments/bas/Analysis/HDDM/best_worst_within_allGS118.pdf", width = 8, height = 4)#, units = 'cm', res = 200, compression = 'lzw'
multiplot( pVDBW,pASVBW, cols=2)
dev.off()

pdf 
  2

testing for residual OV effects

BWRTmod1allin <- lmer(logRT~ cgoalvAvRemBid+cavGoalValue+ cAvBid+isChooseBestTrial+(isChooseBestTrial|subj_idx), data = a118,
                                     REML=FALSE)

sv1_max <- svd(getME(BWRTmod1allin, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(BWRTmod1allin,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("log RT"), pred.labels= c("(Intercept)", "Value Difference", "Overall Goal Value", "Overall Reward Value", "Best - Worst Condition"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.70	1.52 – 1.88	18.77	15.00	<0.001
Value Difference	-0.61	-0.71 – -0.51	-12.05	1630.00	<0.001
Overall Goal Value	-0.21	-0.30 – -0.11	-4.20	1489.00	<0.001
Overall Reward Value	-0.01	-0.11 – 0.09	-0.20	1635.00	0.838
Best - Worst Condition	0.05	-0.10 – 0.20	0.69	15.00	0.499
Random Effects
σ²	0.21
τ₀₀ _{subj_idx}	0.10
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.07
ρ₀₁ _{subj_idx}	-0.30
ICC	0.36
N _{subj_idx}	14
Observations	1644
Marginal R² / Conditional R²	0.066 / 0.404

model comparison

Table S1

BWRTmodbase <- lmer(logRT~ cgoalvAvRemBid+isChooseBestTrial+(isChooseBestTrial|subj_idx), data = a118,
                                   REML=FALSE)

sv1_max <- svd(getME(BWRTmodbase, "Tlist")[[1]]) 
variance_explained <- round(sv1_max$d^2/sum(sv1_max$d^2)*100, 1)  

tab_model(BWRTmodbase,p.val= "kr", show.df = T ,show.stat = TRUE, string.stat="t", dv.labels = c("log RT"), pred.labels= c("(Intercept)", "Value Difference", "Best - Worst Condition"), col.order = c("est", "se","ci",  "stat","df", "p"))

	log RT
Predictors	Estimates	CI	t	df	p
(Intercept)	1.70	1.52 – 1.88	18.73	15.00	<0.001
Value Difference	-0.60	-0.70 – -0.50	-11.78	1628.00	<0.001
Best - Worst Condition	0.08	-0.07 – 0.22	1.08	15.00	0.299
Random Effects
σ²	0.22
τ₀₀ _{subj_idx}	0.11
τ₁₁ _{subj_idx.isChooseBestTrial2-1}	0.06
ρ₀₁ _{subj_idx}	-0.34
ICC	0.36
N _{subj_idx}	14
Observations	1644
Marginal R² / Conditional R²	0.059 / 0.395

xxx <- anova(BWRTmod0, BWRTmod1b, BWRTmod1g, BWRTmodbase)
row.names(xxx) <- c("VD + C<sub>b-w</sub> (baseline)", "baseline + OV<sub>reward</sub>", "baseline + OV<sub>goal</sub>", "baseline + C<sub>b-w</sub>: OV<sub>reward</sub>")
attributes(xxx)$heading <- NULL
xxx2 <- xxx#[, c(1,2, 6, 7, 8)]
xxx2$AIC <- round(xxx2$AIC)
xxx2$BIC <- round(xxx2$BIC)
xxx2$Chisq <- round(xxx2$Chisq,2)
xxx2[,8] <- round(xxx2[,8],3)
xxx2$logLik <- xxx2$AIC - c(NaN, xxx2$AIC[1:length(xxx2$AIC)-1]) 
allR2s <- rsquared(list(BWRTmodbase,BWRTmod0, BWRTmod1b, BWRTmod1g))
xxx2[,1] <- round(allR2s$Conditional,2)
xxx2 <- xxx2[, c(1:4, 6, 8)]
colnames(xxx2) <- c("R<sup>2</sup>", "AIC", "BIC", "dAIC", "Chi<sup>2</sup>", "p")
htmlTable(xxx2)

	R²	AIC	BIC	dAIC	Chi²	p
VD + C_b-w (baseline)	0.33	2255	2293
baseline + OV_reward	0.33	2257	2300	2	0.17	0.676
baseline + OV_goal	0.35	2239	2282	-18	17.77	0
baseline + C_b-w: OV_reward	0.35	2241	2289	2	0.04	0.837

Goal congruency dominates reward value in accounting for behavioral and neural correlates of value-based decision-making

Romy Froemer

2019-08-28

Study 1: best vs worst within subject

accuracy

rt

Liking

LCA simulations:

Goal value LCA

Reward value LCA

Study 2: best vs worst within fMRI

accuracy

RT

Liking

fMRI analyses

1) Bartra ROI

make covariance table for predictors:

Testing whether reward or goal values (or both) explain value network (BOLD) activity

2) Follow up analyses

vmPFC and ventral Striatum separately

Interaction model

Gradient analysis Striatum

Estimates by Regions separately

3) Set Appraisal

Value network

Constituting regions of the value network

Follow up Striatal gradient for set appraisal

4) Supplemental Analyses

PCC

rACC vs mOFC

Study S0: choose best only

accuracy

RT

Study S1: choose worst only

accuracy

rt

Study S2: best vs worst within, incentivized

accuracy

rt